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Abstract 

Evaluation of the eigenvectors of symmetric tridiagonal matrices is one of the most basic 
tasks in numerical linear algebra. It is a widely known fact that, in the case of well separated 
eigenvalues, the eigenvectors can be evaluated with high relative accuracy. Nevertheless, in 
general, each coordinate of the eigenvector is evaluated with only high absolute accuracy. In 
particular, those coordinates whose magnitude is below the machine precision are not expected 
to be evaluated to any correct digit at all. 

In this paper, we address the problem of evaluating small (e.g. 10 -50 ) coordinates of the 
eigenvalues of certain symmetric tridiagonal matrices with high relative accuracy. We propose a 
numerical algorithm to solve this problem, and carry out error analysis. Also, we discuss some 
applications in which this task is necessary. 

Our results are illustrated via several numerical experiments. 

Keywords: symmetric tridiagonal matrices, eigenvectors, small elements, high accuracy, recurrence 
relations 

Math subject classification: 65G99, 65F15, 65Q30 

1 Introduction 

The evaluation of eigenvectors of symmetric tridiagonal matrices is one of the most basic tasks in 
numerical linear algebra (see, for example, such classical texts as [3], [I], [5], [6], [8], [9], [16], [15] . 
[2TJ]). Several algorithms to perform this task have been developed; these include Power and Inverse 
Power methods, Jacobi Rotations, QR and QL algorithms, to mention just a few. Many of these 
algorithms have become standard and widely known tools. 

In the case when the eigenvalues of the matrix in question are well separated, most of these 
algorithms will evaluate the corresponding eigenvectors to a high relative accuracy. More specifically, 
suppose that n > is an integer, that v £ R n is the vector to be evaluated, and v £ M. n is its numerical 
approximation, produced by one of the standard algorithms. Then, 

^ - e, (1) 



where || • || denotes the Euclidean norm, and £ is the machine precision (e.g. e ~ 10~ 16 for double 
precision calculations) . 

However, a closer look at ([T]) reveals that it only guarantees that the coordinates of v be evaluated 
to high absolute accuracy. This is due to the following trivial observation. Suppose that we add 
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e ■ || u|| to the first coordinate v\ of v. Then, the perturbed v will not violate ([T]). On the other hand, 
the relative accuracy of i)i can be as large as 



\vi + e 



H 

M 



(2) 



In particular, if, say, ||w|| = 1 and \v±\ < s, then v\ is not guaranteed to approximate v\ to any 
correct digit at all! 

Sometimes the poor relative accuracy of "small" coordinates is of no concern; for example, 
this is usually the case when v is only used to project other vectors onto it. Nevertheless, in several 
prominent problems, small coordinates of the eigenvector are required to be evaluated to high relative 
accuracy. These problems include the evaluation of Bessel functions (see Sections |3.2M3.4.3"ll6.ip . and 
the evaluation of some quantities associated with prolate spheroidal wave functions (see Sections 13. 31 
16.21 and also [H], [TS]), among others. 

In this paper, we propose an algorithm for the evaluation of the coordinates of eigenvectors of cer- 
tain symmetric tridiagonal matrices, to high relative accuracy. As the running example, we consider 
the matrices whose non-zero off-diagonal elements are constant (and equal to one) , and the diagonal 
elements constitute a monotonically increasing sequence (see Definition [lj . The connection of such 
matrices to Bessel functions and prolate spheroidal wave functions is discussed in Sections l3.4.3[[^^l 
respectively. Also, we carry out detailed error analysis of our algorithm. While our algorithm can 
be viewed as a modification of already existing (and well known) algorithms, such error analysis, 
perhaps surprisingly, appears to be new. In addition, we conduct several numerical experiments, to 
both illustrate the performance of our algorithm, and to compare it to some classical algorithms. 

The following is the principal analytical result of this paper. 

Theorem 1. Suppose that a > 1 and c > 1 are real numbers, and that n > is an integer. Suppose 
also that A\ , . . . , A n are defined via the formula 

A j = 2+(Q , (3) 

for every j — 1, . . . , n, and that the n by n symmetric tridiagonal matrix A is defined via the formula 

\ 



A 



(A x 1 
1 A 2 1 
1 A 3 



V 



l 



l 

A n J 



(4) 



Suppose furthermore that A is an eigenvalue of A, and that 

A > 2 + A k , 



(5) 



for some integer 1 < k < n. Suppose, in addition, that X — (X\, . . . , X n ) € W l is a A— eigenvector 
of A. If the entries of A are defined to the relative precision e > 0, then the first k coordinates 
X%, . . . , Xk of X are defined to the relative precision R(c, a), where 



R(c,a)<e-o(c 4a ^ a + 2 ^ , 



(6) 
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Remark 1. We observe that, according to (JB]), the relative precision of X\, . . . , X% does not depend 
on their order of magnitude. In other words, even if, say, X\ is significantly smaller than e, it is 
still defined to fairly high relative precision. 

Remark 2. The definition of the entries of the matrix A in Theorem]^ is motivated by particular 
applications (see Section^. On the other hand, Theorem^ and Remark^ generalize to a much 
wider class of matrices; these include, for example, perturbations of A, defined via Q , as well as 
matrices, whose diagonal entries are of a more general form than ^ . While such generalizations 
are straightforward (and are based, in part, on the results of Section \4-£\ ), the analysis is somewhat 
long, and will be published at a later date. 

The proof of Theorem [T] is constructive and somewhat technical (see Section 21 in particluar, 
Section B~3l . The resulting numerical algorithm for the evaluation of the eigenvector X is described 
in Sections 15.21 15.31 15.41 and Section [531 contains its error analysis. 

The paper is organized as follows. Section [2] contains a brief informal overview of the principle 
algorithm of this paper. In Section [31 we summarize a number of well known mathematical and nu- 
merical facts to be used in the rest of this paper. In Section [31 we introduce the necessary analytical 
apparatus and carry out the analysis. In Section [SJ we introduce the principal numerical algorithm 
of this paper, and perform its error analysis; we also describe several other related algorithms. In 
Section [6l we discuss some applications of our algorithm to other computational problems. In Sec- 
tion [71 we illustrate the performance of our algorithm via several numerical examples, and compare 
it to some related classical algorithms. 

2 Overview 

In this section, we overview the intuition behind the principal algorithm of this paper and the 
corresponding error analysis. 

Suppose that A is a symmetric tridiagonal matrix, whose non-zero off-diagonal elements are 
constant (and equal to one), and the diagonal elements constitute a monotonically increasing se- 
quence A\ < A2 < ... (see Definition [1]) . Then, the coordinates of any eigenvector of A satisfy 
a certain three-term linear recurrence relation with non-constant coefficients, that depend on the 
corresponding eigenvalue A (see Theorem [7j. More specifically, 



for every k = 2, . . . , n — 1, where x — (x\, . . . , x n ) € R™ is the eigenvector corresponding to the 
eigenvalue A. 

It turns out that the qualitative properties of the recurrence relation (J7j depend on whether 
A — Ak is greater than 2, less than —2, or between —2 and 2. Both our algorithm and the subsequent 
error analysis are based on the following three fairly simple observations. 

Observation 1 ("growth"). Suppose that B > 2 and x%, X2,X3 are real numbers, and that 



If < x\ < X2, then the evaluation of 23 from x\, xi via (J8j) is stable (accurate); moreover, X3 > X2- 
On the other hand, the evaluation of £1 from 22,23 via ([5]) is unstable (inaccurate), since, loosely 
speaking, we attempt to evaluate a "small" number as a difference of two bigger positive numbers. 

Remark 3. This intuition is generalized and developed in Theorem \ll[ Corollary^ Theorem Xl'A and 
serves a basis for the step of the principal algorithm, described in Section \5.3.1\ (see also Observation 
1 in Section \5.5\) . 



Xk+i - (A - A k ) ■ x k + x k -i = 0, 



(7) 



X3 — B ■ xi + £1 = 0. 
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Observation 2 ("decay"). Suppose now that B < —2 and 2/1,2/2,2/3 are real numbers, and 
that 

2/3-5-2/2 + 2/1=0. (9) 

If 2/3 and 2/2 have opposite signs, and I2/2I > I2/3I, then the evaluation of 2/1 from 2/2,2/3 via © is 
stable (accurate); moreover, 2/1 and 2/2 have opposite signs, and |j/i| > 1 2/2 1 - On the other hand, 
the evaluation of 2/3 from 2/1,2/2 © is unstable (inaccurate), since, loosely speaking, we attempt to 
obtain a small number as a sum of two numbers of opposite signs of larger magnitude. 

Remark 4. This intuition is generalized and developed in Theorem \13l and serves a basis for the 
step of the principal algorithm, described in Section \5.3.4\ (see also Observation 3 in Section ] 5. 5]) . 

Observation 3 ("oscillations"). Consider the following example. Suppose that a > is a real 
number, and that the real numbers xq, Xx,x%, ■ ■ ■ , are defined via the formula 

Xk = cos(fc • a), (10) 

for every k = 0, 1, . . . ,. We recall the trigonometric identity 

cos ((k + 1) • a) + cos ((fc — 1) • a) = 2 ■ cos(a) • cos (fc • a) , (11) 

that holds for all real a, k, and substitute (ITU)) into (fTT|) to conclude that 

x k+ i = 2 ■ cos(a) ■ x k - x k -i, (12) 

for every k = 1,2,.... Obviously, the sequence xq, X\, . . . contains elements of order one as well 
as relatively small elements - the latter ones are obtained as a difference of two larger elements, 
potentially resulting in a loss of accuracy, and complicating the error analysis of the recurrence 
relation (fl2)) . 

Consider, however, the complex numbers zq, z\, Zi, • ■ ■ , defined via the formula 

z k = e lka , (13) 

for every k = 0, 1, We combine (fl"U|) and (|13l) to conclude that 

x k = $t(z k ), (14) 

for every k = 0,1,2, Moreover, 

z k +i = e ta ■ z k , (15) 

for every k — 0, 1, ... . The stability of the recurrence relation (fTB"]) is obvious (if its elements are 
stored in polar coordinates, each subsequent element is obtained from its predecessor via rotating 
the latter by the angle a; also, the absolute values of zq, zi, ■ ■ ■ are all the same, in sharp contrast 
to those of xq, x%, . . . ). Moreover, even at first glance, (fT5")) seems to be easier to analyze than (|12l) . 

Remark 5. This intuition is generalized and developed in Theorems \1J\ \15\\1(K\1T\ Corollaries^ 
Theorems \ 23\ \2J\ Corollaries [7] and serves a basis for the steps of the principal algorithm, 
described in Sections 1 5. 3. 3\ Y5. 3. 61 (see also Observations 5-10 in Section \5.5\ ). 

3 Mathematical and Numerical Preliminaries 

In this section, we introduce notation and summarize several facts to be used in the rest of the 
paper. 
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3.1 Real Symmetric Matrices 



In this subsection, we summarize several well known facts about certain symmetric tridiagonal 
matrices. 

In the following theorem, we describe a well known characterization of the eigenvalues of a 
symmetric matrix. 

Theorem 2 (Courant-Fisher Theorem). Suppose that n > is a positive integer, and that A is an 
n by n real symmetric matrix. Suppose also that 



o-i > cr 2 > ■ ■ ■ > a n 

are the eigenvalues of A, and 1 < k < n is an integer. Then, 

cr k = max{min{:E T Az; : x £ U, \\x\\ = l} : U C M", dim(U) = k] , 
for every k = 1, . . . , n. Also, 



(16) 



(17) 



0~k = 

min {max {x T Ax : x G V, \\x\\ = l} : V C R", dim(V) =n — k + l}, (18) 

for every k — 1, . . . , n. 

The following theorem is an immediate consequence of Theorem [5] 

Theorem 3. Suppose that n > is an integer, and A,B are symmetric positive definite n by n 
matrices. Suppose also that 

Oi\ > a 2 > • • ' > Ot n (19) 

are the eigenvalues of A, and that 

Pi > > ■ ■ ■ > Pn (20) 
are the eigenvalues of B. Suppose furthermore that A — B is positive definite. Then, 

Pk < a k , (21) 

for every integer k = 1, . . . , n. 

In the following theorem, we describe the eigenvalues and eigenvectors of a very simple symmetric 
tridiagonal matrix. 

Theorem 4. Suppose that n > 1 is a positive integer, and B is the n by n symmetric tridiagonal 
matrix, defined via the formula 



(2 1 
1 2 1 

1 2 1 



B = 



\ 



1 2 1 
1 2/ 



(22) 



In other words, all the diagonal entries of B are equal to 2, and all the super- and subdiagonal entries 
of B are equal to 1. Suppose also that 

<?i > c 2 > • • • > cr n -i > o- n (23) 

are the eigenvalues of B, and Vi,...,v n €E K™ are the corresponding unit length eigenvectors with 
positive first coordinate. Then, 

a k = 2- (cos (-^r) + l) , (24) 



n + 1 



for each k = 1, . . . , n. Also, 



1 / . / %k \ ( 27rfc \ / rnrk \\ 

Vk = V^' [ sm l^TT J ' sm {— ) ' ' • ■ ' sm U+I J J ' (25) 

/or each k — 1, . . . , n. 

Remark 6. The spectral gap of the matrix B , defined via (|22p . is 

, I _„ 2 = 4 . si »^._^). si »(5._i_) =0 ( ; L), (26) 

dite to (|24[) above. Also, the condition number of B is 

<B) = fi = 1 + ^(^+1)) = 2 
v ; o-„ 1 - cos (7r/(n + 1)) v ' v ' 



3.2 Bessel Functions 

In this section, we describe some well known properties of Bessel functions. All of these properties 
can be found, for example, in [T], [7J. 

Suppose that n > is a non-negative integer. The Bessel function of the first kind J n : C — > C 
is defined via the formula 

J »W = E r? t v • (o) > ( 28 ) 

' m • (m + «)! V2/ 

m=0 v ' 

for all complex z. Also, the function J_ n : C — > C is defined via the formula 

J_„(z) = (-1)" • J n (z), (29) 

for all complex z. 

The Bessel functions Jo, J±i, J±2, ■ ■ ■ satisfy the three-term recurrence relation 

z ■ J n -i{z) + z ■ J n+ x(z) = 2n ■ J n (z), (30) 
for any complex z and every integer n. In addition, 

oo 

E J n(*) = l> (31) 
n— — oo 

for all real x. In the following theorem, we rewrite (|30[) . (|31l) in the matrix form. 
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Theorem 5. Suppose that x > is a real number, and that the entries of the infinite tridiagonal 
symmetric matrix A(x) are defined via the formulae 

A n!n ^x{x) = A n , n+1 (x) = 1, (32) 

for all integer n, and 

2n 

A nn (x) = , (33) 
x 

for every integer n. Suppose also that the coordinates of the infinite vector v(x) are defined via the 
formula 

v n (x) = J n (x), (34) 
for every integer n. Then, v is a unit vector (in the I 2 — sense) , and, moreover, 

A(x) ■ v(x) = 0, (35) 
where stands for the infinite dimensional zero vector. 

3.3 Prolate Spheroidal Wave Functions 

In this section, we summarize several well known facts about prolate spheroidal wave functions. 
Unless stated otherwise, all these facts can be found in [21], [17], [11], [18], [TO] . 

Suppose that c > is a real number. The prolate spheroidal wave functions (PSWFs) corre- 

(c) (c) 

sponding to the band limit c are the unit- norm eigenfunctions ipQ ,ip\ , ■ . ■ of the integral operator 
F c : L 2 [— 1, 1] — > [—1, 1], defined via the formula 

F c [<p](x) = ^ <p{t) ■ e lcxt dt. (36) 

The numerical algorithm for the evaluation of PSWFs, described in [21], is based on the following 
facts. Suppose that Po,P\,... are the Legendre polynomials (see, for example, [T], [?])■ For every 
real c > and every integer n > 0, the prolate spheroidal wave function ipn ^ can be expanded into 
the series of Legendre polynomials 

oo 

i> n (x) = /3i"' c) • Pk (x) ■ VkTTj2, (37) 

fe=0 

for all —1 < x < 1, where P^ ,c \ /3[ n ' c \ . . . are defined via the formula 

Pt ,C) = J ■ p k(x) ■ Vk + T/2 dx, (38) 

for every fc = 0, 1, 2, Moreover, 

(/^ C) ) 2 + (M n,c) ) 2 + (4 n,c) ) 2 + • • ■ = L ( 39 ) 

Suppose also that the non-zero entries of the infinite symmetric matrix A^ are defined via the 
formulae 



k(k | n | 2fc(fc + i)-i , 

h,k~^ + L)+ ( 2 fc + 3)(2fc-l) 



c 



4t +2 = 41 h = (fc + 2)(fc + 1) • c\ (40) 

Kh+2 h+2 ' h (2fc + 3)v/(2fc + l)(2fc + 5) 
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for every k = 0, 1, 2, 

The matrix is not tridiagonal; however, it naturally splits into two symmetric tridiagonal 
infinite matrices A c,even and A c ' odd , defined, respectively, via 



I ' A {c) A [c) 

1 -^0,0 ^0.2 



A c 



.(c) 
1 2,0 



,(c) Ac) 
A 2,2 A 2A 



A 



(c) 



(c) 



4 



and 



.odd 



(a\ 

A 



V 



) 

1.1 

(c) 



Ac) 

A l,3 

4(C) 

^3.3 

43 



.(c) 
^3,5 

.(c) 
^5,5 



.4 



(c) 
5.7 



Suppose that Xo^ < Xjf^ < ■ ■ • arc the eigenvalues of v4 c 
eigenvalues of A c ' odd . If n is even, then 



\ 
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and Xi ] < X ( 3 C) < 



(41) 



.4 <: 



c) 



and = for every odd fc. If n is odd, then 

T 



A 



^■(^ c) ,^ c) ,...) =xl c) -(M"' c) ,^ n ' c) ,-. 



(42) 



are the 



(43) 



(44) 



and /^"'^ = for every even k. The eigenvectors in (14*31 . (l4"4"]l have unit length, due to ([55]l. 

(c) 

The algorithm for the evaluation of starts with computing the corresponding eigenvector of 
A c ' even , if n is even, or of A c ' c , if n of odd. In practice, these matrices are replaced with their 
N x N upper left square submatrices, where N > is a sufficiently large integer (see e.g. [3T] for 
further details) . The coordinates of this eigenvector are the coefficients of the expansion of t/4 into 
the Legendre series (|3"T)l (see also (|43p . (|44j) ) . Once this coefficients are precomputed, we evaluate 
ipn^{x) for any given — 1 < x < 1 via evaluating the sum 



2JV-1 



(45) 



k=0 



in O(N) operations (see e.g. [21], [14], [15] for further details). 
3.4 Numerical Tools 

In this subsection, we summarize several numerical techniques to be used in this paper. 



3.4.1 Power and Inverse Power Methods 

The methods described in this subsection are widely known and can be found, for example, in [3J, 
[B]. Suppose that A is an n x n real symmetric matrix, whose eigenvalues satisfy 

Wi\ > \a 2 \ > \a 3 \ > •■■ > K|. (46) 
The Power Method evaluates o\ and the corresponding unit eigenvector in the following way. 



8 



• Set vo to be a random vector in 1" such that ||uo|| = y/vgVo = 1. 

• Set j = 1 and 770 = 0. 

• Compute Vj = Avj-\. 

• Set rjj = vJ_ l Vj. 

• Set Vj = Vj/\\vj\\. 

• If — 77 j 1 1 is "sufficiently small", stop. 

• Otherwise, set j = j + 1 and repeat the iteration. 

The output value rjj approximates o\ , and Vj approximates a unit eigenvector corresponding to o\ . 
The cost of each iteration is dominated by the cost of evaluating Avj-%. The rate of convergence of 
the algorithm is linear, and equals and the error after j iterations is of order (|o'2|/|ci|)' 7 - 

Remark 7. In this paper, we use the modification of this algorithm, in which i and r)j are evaluated 
via the formulae 



i = argmax{|w :) _i(fc)| : k = !,...,«}, rjj 



(47) 



The Inverse Power Method evaluates the eigenvalue of A and a corresponding unit eigenvector, 
given an initial approximation a of Ofc that satisfies the inequality 

|cr — crfe j < max{|er — a -\ : j ^ k} . (48) 

Conceptually, the Inverse Power Method is an application of the Power Method on the matrix 
B = (A — a I) 1 . In practice, B does not have to be evaluated explicitly, and it suffices to be able 
to solve the linear system of equations 

(-4 - al) i)j = Vj-i, (49) 
for the unknown Vj , on each iteration of the algorithm. 

Remark 8. // the matrix A is tridiagonal, the system (|49|) can be solved in 0(n) operations, for 
example, by means of Gaussian elimination or QR decomposition (see e.g. J3]/ ; JSjj, f20f). 

3.4.2 Sturm Sequence 

The following theorem can be found, for example, in [3D] (see also [2]). It provides the basis for an 
algorithm for evaluating the fcth smallest eigenvalue of a symmetric tridiagonal matrix. 



Theorem 6 (Sturm sequence). Suppose that 



C = 



( a\ 
b 2 





b 2 
a-2 





h 





b n 







6, 
a n J 



(50) 



is a symmetric tridiagonal matrix such that none of bi, ■ ■ ■ ,b n is zero. Then, its n eigenvalues satisfy 

o-i(C) < ••• <<7„(C). (51) 
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Suppose also that Ck is the k x k leading principal submatrix of C , for every integer k — 1, . . . , n. 
We define the polynomials p-i,po, . . . ,p n via the formulae 

p_i(aO = 0, Po(x) = l (52) 

and 

p k (x) = det(C k - xlk) , (53) 

for every k = 2, ...,n. In other words, pk is the characteristic polynomials of Ck, for every k = 
1, . . . , n. Then, 

p k (x) = (a k - x) p k -i{x) - b\pk-i{x), (54) 

for every integer k = 1,2, ... ,n. Suppose furthermore, that, for any real number a, the integer A(o~) 
is defined to be the number of agreements of sign of consecutive elements of the sequence 

Po(o-),Pi(o-),...,p„(a), (55) 

where the sign of Pk{o~) is taken to be opposite to the sign of pk-i{o~) ifpk{o~) is zero. In other words, 

n 

A ( a ) = ^2^ee(p k -i(a),p k (a)) , (56) 
fe=l 

where, for any pair of real numbers a, b, the integer agree(a, b) is defined via 

{1 if ab > 0, 
i ./« = », MO, 
if o = 0, 
if ab < 0. 

Then, the number of eigenvalues of C that are strictly larger than a is precisely A(a). 

Corollary 1 (Sturm bisection). The eigenvalue o~k(C) of (150[) can be found by means of bisection, 
each iteration of which costs 0{n) operations. 

Proof. We initialize the bisection by choosing xq < crfc(C) < yo- Then we set j ' = and iterate as 
follows. 

• Set Zj = (xj + Vj)/2. 

• If yj — Xj is small enough, stop and return Zj . 

• Compute Aj — A{zj) using (f54"| and (|55|) . 

• If Aj > k, set Xj-\-i = Zj and j/j+i = yj. 

• If Aj < k, set Xj+i = Xj and j/j+i = Zj. 

• Increase j by one and go to the hrst step. 

In the end |<7fc(C) — Zj\ is at most yj — Xj. The cost of the algorithm is due to (f5"4")) and the definition 
of A(a). M 
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3.4.3 Evaluation of Bessel Functions 



The following numerical algorithm for the evaluation of Bessel functions (see Section [Q]) at a given 
point is based on Theorem [5] (see e.g pQ). 

Suppose that x > is a real number, and that n > is an integer. The following numerical 
algorithm evaluates Jq(x), J±\, . . . , J± n {x). 

Algorithm: evaluation of Jq(x), J±i, • ■ • , J± n {x). 

• select integer N > n (such TV is also greater than x). 

• set Jjv = 1 and Jn+i = 0. 

• evaluate J/v-i, Jjv-2, ■ ■ ■ , Ji iteratively via the recurrence relation (|30j) . in the direction of 
decreasing indices. In other words, 

2k 



Jk- 



J k {x) - J k +i(x), 



(58) 



for every k = N, . . . , 2. 
evaluate J_i from J\ via ((29)) . 
evaluate Jo from J\, J_i via (|5P|) . 
evaluate d via the formula 



d = 



N 



P 



(59) 



fc=i 



• return Jo/d, J\/d, . . . , J n /d. 

Discussion of such properties of this algorithms as the accuracy to which Jo/d, Ji/d, J n /d 
approximate Jq{x), Ji(x), . . . , J n (x), or the choice of N, are beyond the scope of this paper (see, 
however, Section [7]). Nevertheless, we make the following remark about this algorithm. 

Remark 9. Suppose that x > is a real number, that < n < N are integers, and that the real 
numbers Jo,..., J/v o-rid d are those of the algorithm above (see (|58|) . (|59[) ). Suppose also that A 
is the symmetric tridiagonal 2N + 1 by 2N + 1 matrix, whose non-zero entries are defined via the 
formulae 



A 



.4 



1, 



for every j = 1, . . . , 2N , and 



— Aj — 2 + — , 



for every j — 1, . . . , 2N, 2N + 1. Then, in exact arithmetics, the vector X — {X\, 
M? N+1 , defined via 

X = —■ ^Jjv, • • • , Jl, Jo, — Jti • • • , (— 1)^ ■ Jv) , 
is the unit-length eigenvector of A, corresponding to the eigenvalue X, defined via 

X = 2 + 2 ^ N + 1 \ 



(60) 
(61) 

(62) 
(63) 



In particular, Xi > 0. 
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4 Analytical Apparatus 



The purpose of this section is to provide the analytical apparatus to be used in the rest of the paper. 



4.1 Properties of Certain Tridiagonal Matrices 

In this subsection, we study some general properties of the matrices of the following type. 

Definition 1. Suppose that n > 1 is a positive integer, and < fx < fi < . . . is an increasing 
sequence of positive real numbers. We define the symmetric tridiagonal n by n matrix A = A(f) via 
the formula 



Mi 
i 



V 



i 

1 



1 

A 3 



1 



A n -i 
1 



\ 



1 

An/ 



(64) 



where, for every integer k > 0, the real number A k is defined via the formula 

A k = 2 + f k . 



(65) 



In the following theorem, we describe some properties of the eigenvectors of the matrix A of 
Definition [TJ 

Theorem 7. Suppose that n > 1 is an integer, < f\ < fi < . . . is an increasing sequence of 
positive real numbers, and the symmetric tridiagonal n by n matrix A is that of Definition^ Suppose 
also that the real number X is an eigenvalue of A, and x = (xi,...,x n ) £ R" is an eigenvector 
corresponding to X. Then, 

x 2 = (A - Ai) ■ x\. (66) 

Also, 

x 3 = (X - A 2 ) ■ x 2 - Xl, 
Xk+i = (A - Ah) ■ x k - x k -i, 



x n = (A — A„_i) • %n-i — Xn-2- (67) 

Finally, 

x n -i = (A - A n ) ■ x n . (68) 

In particular, both x\ and x n differ from zero, and all eigenvalues of A are simple. 

Proof. The identities (|6"6")l , (|6"T|) . (|68p follow immediately from Definition [T] and the fact that 

Ax = Xx. (69) 

To prove that the eigenvalues are simple, we observe that the coordinates x 2 , ■ ■ ■ ,x n of any eigen- 
vector corresponding to A are completely determined by x\ via ([66]) . (|67| (or, alternatively, via 
(|68p . (|67p ). Obviously, neither x\ nor x n can be equal to zero, for otherwise x would be the zero 
vector. ■ 
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The following theorem follows directly from Theorem [7J and Theorem [5J in Section 13.4.21 

Theorem 8. Suppose that n > 1 is an integer, < /i < fa < . . . is an increasing sequence 
of positive real numbers, and the symmetric tridiagonal n by n matrix A is that of Definition [7J 
Suppose also that 1 < k < n is an integer, that A& is the kth smallest eigenvalue of A, and that 
x = (xi, . . . , x n ) € K™ is an eigenvector corresponding to A&. Then, 

n-l 

^ agree (xj, = - 1, (70) 

where agree is defined via (|57|) in Theorem^ 

Proof. We define the n by n matrix C via the formula 

C = \ k -I-A. (71) 



Suppose that, for any real t, the sequence Po(t), . . . ,p n (t) is the Sturm sequence, defined via (|55|) 
in Theorem [S] Suppose also that x = (x%, . . . , x n ) is the eigenvector of A corresponding to and 
that 

xi = 1. (72) 

We combine {72J with JBB}, ([57} of Theorem [7J and {52}, ((SJl) of Theorem [5] to conclude that 

Po(0) = a;i ) pi(0) = x 2 , Pn-i(0) = x n . (73) 

Moreover, due to the combination of (|r?5|) and flM} , 

P«(0) = 0. (74) 

We observe that the matrix C, defined via (J7TJ), has k — 1 positive eigenvalues, and combine this 
observation with {55]), (J37J) of Theorem [5] and ([75 ]) . ([73 ]) to obtain (JTUJ). ■ 

Corollary 2. Suppose that n > 1 is an integer, < fx < fa < ■ ■ ■ is an increasing sequence of 
positive real numbers, and the symmetric tridiagonal n by n matrix A is that of Definition^ Suppose 
also that Ai and A n are, respectively, the minimal and maximal eigenvalues of A. Then, 

Ai < Ax, (75) 

and 

A„ > A n . (76) 



Proof. Suppose that x — (x\, . . . , x n ) is an eigenvector corresponding to Ai. Due to the combination 
of Theorem [JJ and Theorem [51 

xx ■ x 2 < 0. (77) 

We combine ([77j) with (|66p to obtain ([75)1 . On the other hand, due to the combination of Theorem[7J 
and Theorem [S[ 

X n -i ■ x n > 0. (78) 

We combine (J7SJ with (J55J to obtain ([75]). ■ 
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The following theorem is a direct consequence of Theorems [3] S] in Section 13.11 

Theorem 9. Suppose that n > 1 is an integer, < f% < < . . . is an increasing sequence 

of positive real numbers, and the symmetric tridiagonal n by n matrix A is that of Definition QJ 
Suppose also that 

Ai < A 2 < ■ • • < A» (79) 

are the eigenvalues of A. Then, 

A fe >2- (l-cos(^- ] ] , (80) 



n + 1 

/or every integer k = 1, . . . , n. In particular, A is positive definite. Also, 

A fe > A fe - 2, (81) 

/or every k = 1, . . . , n. 

Proof. The inequalities ([50)1 , (|5Tj) follow from the combination of Theorems [21 HI in Section 13.11 and 
Definition [TJ ■ 

In the following theorem, we provide an upper bound on each eigenvalue of A. 

Theorem 10. Suppose that n > 1 is an integer, < fx < fi < . . . is an increasing sequence 
of positive real numbers, and the symmetric tridiagonal n by n matrix A is that of Definition [JJ 
Suppose also that 1 < k < n is an integer, and that Afe is the kth smallest eigenvalue of A. Then, 

Afe < 2 + A k . (82) 

Proof. Suppose, by contradiction, that 

Afe > 2 + A k . (83) 
Suppose that x = (x%, . . . , x n ) € 1" is an eigenvector of A corresponding to Afe. Suppose also that 

xi > 0. (84) 

It follows from ([83]) that 

Afe - A x > Afe - A 2 > ■ ■ ■ > Afe - A k > 2. (85) 
We combine ([66]), ([67]) in Theorem [7] with (|M]l. ([85]) to conclude that 

x x > 0, x 2 > 0, Xfc + i > 0. (86) 

We combine ([86]) with ([57]) in Theorem [6] in Section l"3.4. 21 to conclude that 

^2 a § ree (a^ , Xj+i ) > fc, (87) 
in contradiction to Theorem [5] ■ 
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4.2 Local Properties of Eigenvectors 

Throughout this subsection, n > is an integer, and < f± < fa < . . . is an increasing sequence 
of positive real numbers. We will be considering the corresponding symmetric tridiagonal n by n 
matrix A = A(f) (see Definition [1] in Section ETTj) . 

In the following theorem, we assert that, under certain conditions, the first element of the eigen- 
vectors of A must be "small" . 

Theorem 11. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] 
in Section \4-l\ Suppose also that X is an eigenvalue of A, and that x — {x\, . . . ,x n ) € R™ is a 
corresponding eigenvector whose first coordinate is positive, i.e. X\ > 0. Suppose, in addition, that 
1 < k < n is an integer, and that 



X>A k + 2. 
Then, 

< xi < x 2 < ■ ■ ■ < Xk < Xk+i- (89) 

Also, 



for every j = 2, . . . , k. In addition, 



*j ^ A " A 3 ■ , . ( A - A 3 



1 ^ ^ - ^. (91) 

X k -i X 2 Xi 



Proof. It follows from d88|) that 

A fe - A x > X k - A 2 > ■ ■ ■ > X k - A k > 2. (92) 

We combine (p7)l in Theorem [7] with (l9"2l to obtain (|89p by induction. Next, suppose that the 
real numbers ri , . . . , are defined via 

r, = ^±i, (93) 

Xj 

for every j = 1, . . . , k. Also, we define the real numbers o\, . . . , o~k via the formula 



for every j = 1, . . . , k. In other words, <7j is the largest root of the quadratic equation 

x 2 -(A- Aj)-x + 1 = 0. (95) 

We observe that 

o-i > ■ ■ ■ > a k > 1, (96) 
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due to {SI]) and EH- Also, 

ri > <Ti > cr 2 > 1, (97) 

due to the combination of (|94[) and (|66p . Suppose now, by induction, that 

r,-_i > (7 j > 1. (98) 

for some 2 < j < fc — 1 . We observe that the roots of the quadratic equation ([55)1 are 1 /<7j < 1 < aj , 
and combine this observation with (|58"|) to obtain 



r?_i - (A - A,-) • r,-_i + 1 > 0. 
We combine (EH with EH) and EH) to obtain 



x 



j+i _ (A - Aj) ■ Xj - Xj-i 



A- A 



< r,-_i. 



3 r. i J 



Also, we combine (|M|). (|58"]h (IT00)) to obtain 

1 



, j — A - A,- - - > A - Aj - — — - — u 3 ^ ^ 

In other words, implies (|101[) . and we combine this observation with (1^71) to obtain 

ri > <7 2 , r 2 > cr 3 , r fc _i > <7 fc . 

Also, due to (fTOOl) . 

ri > r 2 > • • • > r-k-x. 
We combine O, ® 5 ([HE]) . (fTU51) to obtain (jgDJl, (|M|). 
Corollary 3. Under the assumptions of Theorem MR 



1 (A - Aj) • (jj 1 



(99) 
(100) 

(101) 

(102) 
(103) 



Xx 



> 



J=2 



A ~ Aj 



- 1 



(104) 



Remark 10. In J13f the derivation of an upper bound on the first coordinate of an eigenvector 
of a certain matrix is based on analysis, similar to that of Theorem 

Theorem 12. Suppose that the integers 1 < k < n, the real numbers A\, . . . ,A n and X, and the 
vector x = [x\, . . . , x n ) € R™ are those of Theorem \ll\ Suppose also that 2 < j < k is an integer, 
and that £j—i, Sj are real numbers. Suppose furthermore that the real number Xj+x is defined via the 
formula 



Then, 



In particular, 



x j+ x = (A - A 3 ) ■ Xj ■ (1 + Sj) - Xj-x ■ (1 + £j-x) 



Xj + X - Xj + X 



x j+1 



Xj+x - Xj+X 
Xj+X 



< Si 



A — Aj 
A - Aj - 1 



2 + A k - A 3 

1 + A k ~ Aj 



(x-Aj-iy 



(l + A k -AjY 



(105) 
(106) 

(107) 
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Proof. We combine (I105P with ([fJT|) of Theorem [7] to obtain 

x j+x - x j+1 = (A — A,-) • lEj • £j - xj-i ■ £j-x. (108) 
Also, due to the combination of (p7| and of Theorem [TT1 

<^^"T- (109) 



Xj + i A - Aj - 1 
Then, due to {STJ of Theorem [III 

— . (no) 

Xj+i x j+i 

We combine (fID5]l . pTO|) and (fTTU| to obtain (TITO)) . We combine ([TUrJl) with flSSJl to obtain (flDT) . ■ 

Remark 11. While the analysis along the lines the proof of Theorem \ll\ would lead to a stronger 
inequality that (|106|) . £/ie latter is sufficient for the purposes of this paper. 

In the following theorem, we study the behavior of several last elements of an eigenvector of A. 

Theorem 13. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] 
in Section \4-l\ Suppose also that A is an eigenvalue of A, and that x — {x\, . . . , x n ) € R™ is a 
corresponding eigenvector whose last coordinate is positive, i.e. x n < 0. Suppose, in addition, that 
1 < k < n is an integer, and that 

\<A k - 2. (Ill) 

Then, 

< \x n \ < \x n -i\ < ■ ■ ■ < \x k \ < \x k -i\. (112) 

Also, 



Xj Aj — X (X-A^ 1 



Xj+l 

for every j = k, . . . , n — 1. In addition, 



>-hr^ + \ \-^r L ) -i. (H3) 



.!> _£^> ...>E2za>E2zl. (114) 



Proof. We defined the real numbers yi, . . . , y n via the formula 

Vj = (-l) i+1 • Zn+l-j, (H5) 

for every j = 1, . . . , n. Also, we define the real numbers Bi, . . . , B n via the formula 

= -4+Hi (H6) 
for every j = 1, . . . , n. In addition, we define the real number fi via the formula 

H = -A. (117) 
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We combine (flTS) , (|TTC| . (fTT7|) with §7$, (JHEJ) of Theorem [7] to obtain 

Vi = (M - Si) • yi > (118) 

and 

= (M - b j) ■ Vj - Vj-i; (H9) 
for every j = 2, . . . , n — 1. We combine (|116p . (j 11 7|) with (| 1 1 1[) to obtain 

M > Bn+i-fc + 2. (120) 

Following the proof of TheoremQl] (with Xj , Aj , A replaced, respectively, with ?/j , , /i) we use (|118p , 
prgj) . (TTZQ]) to obtain 

< J/l < 2/2 < • • • < Un+l-k < Vn+2-k, (121) 

as well as 



2 

for every j = 2, . . . , n + 1 — k, and 



(122) 



1<^±1^<...<M < ^. (123) 
J/n-fe 2/2 2/1 

Now (HU), dH, (HI follow from the combination of (fTT5|) . (fTTSjl . (fTTTl) . (TT2T1) . (Tmi) and (Tmi) . 

■ 

The following definition will be used to emphasize that the elements of an eigenvector of A exhibit 
different behaviors "in the beginning" , "in the middle" and "in the end" . 

Definition 2. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] 
in Section \4-l\ Suppose also that X is an eigenvalue of A, and that x = {x\, . . . ,x n ) € W 1 is a 
corresponding eigenvector whose first coordinate is positive, i.e. x\ > 0. The eigenvalue X induces 
the division of the set {1, . . . ,n} into three disjoint parts Ig(X), Jo (A), Id{X), defined via the formulae 



/g(A) ={l<j<n 
/o(A) ={1< j<n 
I D {X) ={l<j<n 



A - Aj > 2} , (124) 
-2 > A - Aj > 2} , (125) 
-2>X-Aj}. (126) 



In other words, if 

X-A k >2>X- A k+1 > ■ •• > A - A m > -2 > X - A m+1 (127) 
for some integer k,m, then 

I G (X)={l,...,k}, 
Jo(A) ={k + l,...,m}, 

I D (X) = {m + l,...,n}. (128) 

The sets Ig(X), Jo(A), Id(X) are referred to as the "region of growth", "oscillatory region" and "re- 
gion of decay", respectively. 
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Remark 12. Obviously, depending on X, some of /g(A), io(A), /d(A) may fee empty. For example, 
for the matrix A of Theorem^ in Section \3.1l both Ig{X) and -To (A) are empty for any eigenvalue 
A. 

Remark 13. The terms "region of growth" and "region of decay" are justified by Theorems llll [731 
respectively. Loosely speaking, in the region of growth, the elements of x have the same sign, and 
their absolute values grow as the indices increase. On the other hand, in the region of decay, the 
elements of x have alternating signs, and their absolute values decay as the indices increase. 

In the rest of this subsection, we investigate the behavior of the elements of an eigenvector 
corresponding to A in the oscillatory region Io{X) (see Definition [2]). In the following theorem, we 
partially justify the term "oscillatory region" of Definition [2] 

Theorem 14. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [JJ 
in Section \4-l\ that A is an eigenvalue of A, and that x = (x\, . . . ,x n ) G R™ is a corresponding 
eigenvector. Suppose also that the set Io(X), defined via (|125l) in Definition® is 

Io(X) = {k+l,...,m}, (129) 

for some integer 1 < k < m < n. Suppose furthermore that the real numbers 9k, ■ ■ ■ , 9m— l o^e defined 
via the formula 



A - Aj- 



(130) 



2 

for every j = k, . . . , m — 1, that the two dimensional vectors Vk, ■ ■ ■ , v m are defined via the formula 

I Xj I 131. 



K X 3 + K 

for every integer j — k, . . . ,m, and that the 2 by 2 matrices Qk, ■ • ■ , Q m -i cire defined via the formula 

^=(-1 2 003(0;))' (132) 
for every integer j = k, . . . ,m — 1. Then, 

Vj+x = Qj ■ Vj, (133) 

for every integer j = k, . . . , m — 1. 



Proof. We observe that 



< 9 k < 9 k+ i < ■ ■ ■ < 9 m -i < 7T, (134) 



due to the combination of (|125[) and (|130l) . We also observe that the recurrence relation (|133[) follows 
immediately from the combination of Theorem [7] and (|130[) , (|131l) , (|132l) . ■ 

In the following theorem, we describe some elementary properties of a certain 2 by 2 matrix. 

Theorem 15. Suppose that < 6 < tt is a real number. Then, 



-1 \ / 1 \f 1 e M \ fe w 

2sin(0) \~l e- i6 )\-l 2cos{9))\e ie l) \0 e~ u 



(135) 
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Suppose also that x,y are real numbers. Then, 



K V J V s 

where the complex numbers a, fi are defined via 



,i a 0°) +p ( e i)> (i36) 



-1 \ fx 



PJ 2sin(6») V - 1 e J \y 



(137) 



moreover, 

P = e" m ■ a, (138) 

where a denotes the complex conjugate of a. In particular, 

x = 2R(a), (139) 
y = 23?(e ie • a), (140) 

where 3?(z) denotes the real part of any complex number z. 

Proof. The proof is straightforward, elementary, and is based on the observation that 



1 e w \ 1 % fe~ ld -1 



1 / 2sin((9) V _1 e 



(141) 



In the following theorem, we evaluate the condition number of the 2 by 2 matrix of Theorem II 51 

Theorem 16. Suppose that < 9 < -k is a real number, and that n(Q) is the condition number of 
the matrix 

m = Qe f). (142) 

In other words, k{9) is defined via the formula 



<e) = ^§-, (143) 

CT2 (V) 



where o~i(9),o~2(9) are, respectively, the largest and smallest singular values ofU{9). Then 

Proof. The proof is elementary, straightforward, and is based on the fact that the eigenvalues of the 
matrix U(9) ■ U(9)* are the roots of the quadratic equation 

x 2 - 2x + sin 2 (9) = 0, (145) 

in the unknown x. ■ 

In the following theorem, we introduce a complexification of the coordinates of an eigenvector of 
A in the oscillatory region (see Definitions [1] [5]) . 
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Theorem 17. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] 
in Section \4-l\ that A is an eigenvalue of A, and that x — (x%, . . . , x n ) G R n is a corresponding 
eigenvector. Suppose also that the set IoW, defined via (I125[) in Definition^ is 



Io(X) = {k + l,...,m}, 



(146) 



for some integer 1 < k < m < n. Suppose furthermore that the real numbers 6k, ■ ■ ■ , 6m-l are defined 
via (|130[) of Theorem ] 14\ and that the complex numbers ctk, • • • , ot m -\ and /3k, ■ ■ ■ , fim-i « r e defined 
via the equation 



1 e i6 i\ fa 



x j+1 



in the unknowns ctj, /3j, for j = k, . . . , m — 1. Then, 

Xj — 23? (ctj ) , 



,•+1=23* (,. ..<•" ) 



(147) 



(148) 
(149) 



for every j — k, . . . , m — 1. Moreover, 
a j+ i = 



cos(9j) — cos(9j + i) 



e 3 . P rft+ 9 i+i 



sm 



) ■ sin 



(150) 



for every j = k, . . . , m — 2. 



Proof. The identities (|148[) , (I149p follow from the combination of Theorems [TU [15] To establish 
(|150[) . we proceed as follows. Suppose that k < j < m — 2 is an integer. We combine (|131[) . (|132p . 
(fH3"j) of Theorem [H with (TUT)) to obtain 



x j+1 

X j+2 







1 e ie A fa 



-1 2cos(0j) / 1 
We substitute ([135]) . ([HI]) of Theorem [15] into (|T5T1) to obtain 



(151) 



x j+ A _ ( 1 e ie A /e^ 
/ ~ \e lB > 1 / V e ~ l<9 



1 ( "' \ / a 3 ■ e ie * 
e i9 i 1 ) [fSj ■ e 



(152) 



Next, we combine (fT55j) . ([147]) and ([132]) to obtain 

-l 



a j+1 



1 e ^+i\ * / 1 e ie A / ttj . e ^ 
1 / 1 - e i0 i+i\ ( 1 e ie * 

1 _ e 2i8 j+1 y- e i6j+i 1 J [ei0j 1 

1 /l — e i( - ej+1+e ^ e i8j — e i9j+1 



I _ e 2iG j + 1 I e »0j _ giflj+i i _ gi(9 j+ i+fij) 



■ ■ e 

Qi • e 



-US, 



(153) 
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It follows from (|153[) that 



-,-2iB; 



16a . 3 

ctj ■ e 3 H 



^ g2i#j + i 



. q7-. e -»(9i+9j+i) . ( e iSj+i _ gi(20 j+ i- 



= a,- ■ e 



1 g2z#j + i 



a,- • e l9 ' + 2i ■ 3 



^ — g2«#j+i 

. . e i(8]+8j+i) . ( e <«i+i _ e ^) 

1 — g2i$j + i 



1 - e - 2t9 J+ 1 



where ^(z) denotes the imaginary part of a complex number z. We observe that, for 

l-e" sin(t/2) 



numbers i, /i, 



exp 



l-e 4/l sin(/i/2) 
and use (|155[) to obtain 

1 _ e »(e^e J+1 ) S in((g j - -g j+1 )/2) 



1 — e 2»0j+i 



Due to (jT5eJj) . 



sin(0 i+ i) 
e i(ej+e J+ i) , ^gie j+ i _ e ie^ 

^ g2i#j + i 

sin((^-^ +1 )/2) 



• exp 



(Oj - 3 • 6< i+ i) 



sin(^ +1 ) 



• exp 



7j f 



However, 



*7 _ p jr+l 



cos(^j) — cos(0j + i) 
2 y ' 2 • sin ((0,- + e j+ i)/2) 

Finally, we substitute ([15T)) . ([155)) into (ITM| to obtain ([150)) . 

Corollary 4. Under the assumptions of Theorem \17\ 

cos(9j) — cos(0j+i) 



ckj+1 - ay • e ta i\ < |ay • 

sin(0 i+1 ) ■ sin f «a22l±i 
/or every j = k, . . . , m — 2. 

Corollary 5. Under the assumptions of Theorem \ll\ 

x j + = 4i? i ' ( 2 + cos (%) ' cos ( e j + 2 0)) > 
/or every j = k, . . . , m — 2, where the real numbers Rj > and £j are defined via 



ay = Rj ■ ^ 



for every j = k, . . . , m — 2. 



(154) 
any real 

(155) 
(156) 



(157) 



(158) 



(159) 



(160) 
(161) 
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Proof. We combine (I147P of Theorem H7I and (|138[) of Theorem [15] to obtain 

x] + x 2 j+l = 2 • flay | 2 + |/3, | 2 + cos(^) -{aj-Tj + aj- ft")) 
= 2 ■ (2 • |aj| 2 + 2 • cos(^) ■ (a, • /3~)) 

= 4 • flc^ 2 + cos(0,) • 3? (a 2 • e^)) . (162) 

We substitute (fTUTj) into ([T62]) to obtain (fT60|) . ■ 

In the following theorem, we discuss the "zero-crossing" of certain sequences. 

Theorem 18. Suppose that k > is an integer, and that 

< 0i < 6» 2 < ■ • • < 9 k < it (163) 

are real numbers. Suppose also that e > is a real number, that the sequence yi, ■ ■ ■ ,yk+2 is defined 
via the formulae 

Vi =V2 = 1 (164) 



Vi+A = (0 1 \ Vj 



(165) 



for every j = 1, . . . , k, and that the sequence z\, . . . , Zfc+2 is defined via the formulae 

zi = l, z 2 = l + e (166) 

and 



Zj+ A = o l \ Zj 

Z3+2J \-l 2cos(6»j)y \zj+i 



(167) 



for every j — 1, . . . , k. Suppose furthermore that 

yi ,...,y k+1 >0. (168) 

Then, 

Zi, . . .,z k +i > 0, (169) 
In other words, the sequence z\, . . . , z k + 2 cannot "cross zero" before the sequence 3/1, ... , yu+2 does. 
Proof. It follows from (|168p that 

V2 Vk+i 



> 0. (170) 



Also, we combine ATM]) . (IT65)) . (fT66)l . (fT67j) to obtain 



and the recurrence relations 



yi = l < l + e=^, (171) 



l^±a=2cos(^)--^-, (172) 
V3+1 Vj+i ' 



23 



for every j = 1, . . . , m, and 



Zj+l Zj+l 

for every j = 1, . . . ,m. Suppose now that 1 < I < k, and that 



Zj+2 =2cos(^)-— , (173) 



yj±i < fw. ( i 7 4) 

2/i ^ 

We combine ([T72]) . ([175]) . (jTT4"l) to obtain 

= 2cos(# ; ) < 2cos(6>/) = . (1'5) 

m+i yi+i zi+i zi+i 

In other words, (|174[) implies (1175[) . and we combine this observations with (|171l) to obtain, by 
induction on I, 

y±±i < z j±i i (ire) 

yi zi 

for all I = 1, . . . , k. We combine (TITO]) with ([175]) to obtain (ITM| . ■ 

The following theorem is similar to Theorem [TBI above. 
Theorem 19. Suppose that k > is an integer, and that 

0<6 1 <6 2 <---<6 k <9<iT (177) 
are real numbers. Suppose also that the sequence yi,. ■ ■ ,yk+2 is defined via the formulae 

2/1=2/2 = 1 (178) 

and 



(vs+A = ( 1 )( Vi 
\yj+2j \-i 2cos(^)y Vz/j+i 



(179) 



for every j = 1, . . . , k, and that the sequence t%, . . . , is defined via the formulae 

h=t 2 = l, (180) 

and 



t j+ A = ( o l \ / ^ 

t t+2 y i-i 2cos(0)y u i+1 



(181) 



for every j — 1, . . . , k. In other words, as opposed to (|179|) . the matrix is (|18ip is the same for every 
j . Suppose furthermore that 

t u ...,t k+1 > 0. (182) 

Then, 

yi ,...,y k+1 >0, (183) 
In other words, the sequence yi, . . . , j/fc+2 cannot "cross zero" before the sequence ti,..., t k +2 does. 



24 



Proof. The proof is almost identical to the proof of Theorem [TBI an d is based on the observation 
that 

2-cos(0,-) - - > 2-cos(0) - -, (184) 
p r 

for any j = 1, . . . , k, and any real numbers < r < p. ■ 

In the following theorem, we address the question how long it takes for a certain sequence to 
cross zero. 

Theorem 20. Suppose that k > and I > are integers, and that 

o<e k < e k+1 <■■■< e k+l < ix (185) 

are real numbers. Suppose also that e > 0, and that the sequence x k , ■ ■ ■ ,x k +i+2 is defined via the 
formulae 

x k = 1, x k+ i = 1 + e, (186) 



x j+1 \ = f 1 



(187) 



for every j = fc, . . . ,k + I. Suppose furthermore that 



l-\yO k+ i<^. (188) 
Then, 

x k ,x k+ i,...,x k+ i > 0. (189) 

Proof. Suppose that the sequence t k , ... , t k+ i + 2 is defined via the formulae 

t k = t k+1 = 1 (190) 

and 



t j+1 \ = ( 1 \(t j 

t j+2 J l-l 2cos(e k+l )J Uj+xJ' 



(191) 



for every j — k, . . . , k + I (note that the matrix in (|191[) is the same for every j). It follows from 
(ITM| . (|T5T|) that 



(t k+m \ = (0 1 \"Vl 

^fe+m+iy W 2cos(0 fc+i ); 

for every m — 0, . . . ,1 + 1. We observe that 

1 \ m 1 /l e 



-1 2cos(6»W l- e 



(192) 



(193) 
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for every real 9 and every integer m — 0,1, . . . (see also Theorem ITS]) . We combine (|192l) . (|193p to 
obtain 

t k+m \ 1 ( 1 e iefe +'\ / e m9k + l 



In other words, 



tk+ m +i J 1 + e i9 *+' V e +1 1 / \ e ~ m +l 

e -iS k+i /2 / e imG k + l , e -j(m~l)e fc+1 



2cos(6» fe+; /2) \ e < m+1 '> 9 ^ + e - im6 "+' 

1 (cos ((to - 1/2) ■ 9 k+ i) 

cos(0 fe+/ /2) ^cos((TO + l/2)-0 fc+i ) 



(194) 



cos((TO-l/2)-g fc+; ) 

^fc+m = 7^ 7^7 , (195 

COS [Ok+l/2) 

for every m = 0,...,i + 2. We combine (|195[) with (|188[) to conclude that 

**+!,..., tjk+i >0. (196) 

We combine fT90]t . (|T9T1) . (|T96l) with Theorems [H [H to obtain (fT89|) . ■ 

Corollary 6. Suppose, in addition, that I > 3, and £/ia£ k + l<j<k + l — 2 is an integer. Then, 

xj+i > ? > ^r 1 - (197) 



Proof. We combine (|T551) . ([1871) and (lTiS9f to obtain 

Xj +2 — 2xj + i cos(9j) — Xj > 0, (198) 

and hence 

We combine ([TBTl) . (|T59l and (ITM)) to obtain 

-y < = 2a^j cos(#j_i) — Xj-i, (200) 

and hence 

/ 1 \ 3x 

< Xj ■ ( 2cos(0,_ 1 ) --j<^. (201) 

We combine (fT9"9")) . (|2lTT|) to obtain (fTST)) . ■ 

The following theorem follows from the combination of Theorems QTl [T71 [2D] and Corollary El 

Theorem 21. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [JJ 
in Section \4-l\ that X is an eigenvalue of A, and that x = {x\, . . . ,x n ) G R™ is a corresponding 
eigenvector with positive first coordinate, i.e. x% > 0. Suppose also that 

\-A 1 >--->\-A k >2>\- A k+1 > • • ■ > A — A fe+;+1 > 0, (202) 
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for some integer 1 < k < n and 1 <l <n — k — 1. Suppose furthermore that 

A - A k+l + i / 7T 



2 >cos l^rj- (203) 

Traen, 

< xi < x 2 < • • • < x k < x k +i- (204) 

Also, 

x k ,x k+1 ,...,x k +i > 0. (205) 

//, in addition, I > 3, £/ien 

. 2 ' x k 2 ■ x fc+ i 2 ■ x k+ i-3 , , 

^fc+i > , x fc+ 2 > , x fc+i _ 2 > 5 , (206) 



Xk+i-i > — x — • (207) 



Proof. We obtain (|204[) from the combination of (|202p and Theorem QT] Also, we obtain (|205p from 
the combination of Theorems [T71 |2T)1 Finally, we combine Theorems [T71 |2"01 and Corollary [5] to obtain 

PSD , HMD - ■ 

Theorem 22. Suppose that the integers 1 < k < n and I > 3, the real numbers A\, . . . , A n and A, 
and £/ie vector x — (xi, . . . ,x„) € W 1 are those of Theorem \21\ Suppose also that fc + 1 < j < k + l — 3 
zs an integer, and that Ej-\,Ej are real numbers. Suppose furthermore that the real number Xj+i is 
defined via the formula 

x j+1 = (A - Aj) ■ x 3 • (1 + e 3 ) - Xj_! • (1 + ej _!) . (208) 

T/ien, 



Xj+i - Xj+i 



Xj+i 



Proof. We combine (|105l) with (1571) of Theorem [7J to obtain 

Xj+i ~ Xj+x = (A - A,-) • x 3 • £j - Xj-i ■ Ej-i. (210) 
Also, due to the combination of (|67l) and (I206P of Theorem |2~T1 

< £. (211) 

xj+i 4 

We combine (j2TU|) . ([2TT1) to obtain (|2"0gi) . ■ 

Remark 14. While a more careful analysis would lead to a stronger inequality that (|209[) . the latter 
is sufficient for the purposes of this paper. 



27 



The following theorem is almost identical to Theorem Q3J with the only difference that the 
recurrence relation is reversed. 

Theorem 23. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [JJ 
in Section \4-l\ that A is an eigenvalue of A, and that x — (x\, . . . , x n ) G R n is a corresponding 
eigenvector. Suppose also that the set Jo (A), defined via (If 25)) in Definition^ is 

IoW = {k + l,...,m}, (212) 

for some integer 1 < k < m < n. Suppose furthermore that the real numbers ifk, ■ ■ ■ , tpm—~L ar £ 
defined via the formula 

Lfj = arccos ( ^ J+1 J , (213) 

for every j = fc, . . . , m— 1 , that the two dimensional vectors Wk, ■ ■ ■ , ui m+ i are defined via the formula 

for every integer j = k, . . . ,m + 1, and that the 2 by 2 matrices Pk, . . . , Pm—l are defined via the 
formula 

P > = (-1 2co S fe))' <215) 
for every integer j — k, . . . ,m — 1. Then, 

wj = Pj_! ■ w j+ i, (216) 

for every integer j = k + 1, . . . , m. 

Proof. The proof is essentially identical to that of Theorem [TU ■ 

The following theorem is similar to Theorem 1171 with the only difference that the recurrence 
relation is reversed. 

Theorem 24. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] 
in Section \4-l\ that X is an eigenvalue of A, and that x — (x\, . . . , x n ) £ M" is a corresponding 
eigenvector. Suppose also that the set Io{X), defined via (I125p in Definition^ is 

Io(X) = {k + l,...,m}, (217) 

for some integer 1 < k < m < n. Suppose furthermore that the real numbers tfk, ■ ■ ■ , <Pm—~L are 
defined via (12131) of Theorem \23l and that the complex numbers 7fc, . . . ,7 m -i and <5fc, • • • > <5m— 1 are 
defined via the equation 



1 e*A hj\ = {x j+2 -(-ir-l- 2 \ 



(218) 



in the unknowns jj, Sj, for j = k, . . . , m — 1. Then, 

x j+2 =2M( lj )-(-ir-3- 2 ) (219) 
x j+ i = 23? ( 7j • e 1 ^) • (-I)™"'"- 1 , (220) 
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for every j — k, . . . ,m — 1. Moreover, 
7i-i = 

C0s((/3j) — cos(yjj_i) 



■7j 



/or every j = k + 1 , . . . , m — 1 . 

Proof. The proof is essentially identical to that of Theorem \I7\ and will be omitted. 
Corollary 7. Under the assumptions of Theorem \24\ 

cos(ifj) — COs(<fij-i) 



U'-i-7i-e^ < | 7 .| 



sin(^_ 1 ).sin(^±f^) 



(221) 



(222) 



(223) 



(224) 



/or every j = k + 1 , . . . , m — 1 . 

Corollary 8. Under the assumptions of Theorem \24\ 

x 2 j+1 + x 2 j+2 = 4R 2 • (1 + cos(^) • cosOj + 2&)) , 
for every j = k, . . . , m — 1, where the real numbers Rj > and £j are defined via 

for every j = k, . . . , m — 1. 

Proo/. We combine (|218j) of Theorem |2"H and (|138j) of Theorem [131 to obtain 

^ 2 + i+^ 2 + 2 =2- (| 7j | 2 + |^| 2 + cos(^)- (tj-^ + TJ-^)) 
= 2 ■ (2 ■ | 7 , | 2 + 2 ■ cos(^) • K ( 7i • 

= 4.(| 7j | 2 + cos(^)-^(7 2 -e^))- (225) 

We substitute ((224)) into (|225|) to obtain ([223]) . ■ 

The following theorem is analogous to Theorem [5T] 

Theorem 25. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] 
in Section \4-l\ that A is an eigenvalue of A, and that x = (x±,...,x n ) G R™ is a corresponding 
eigenvector. Suppose also that 



> A - A m - r > • • • > A - A m > -2 > A - A m+ i, 
for some integer 1 < m < n and 1 < r < m — 1. Suppose furthermore that 



Am—r A 



> COS 



2r - 1 / ' 



and that x m > 0. Then, 

< x n ■ (-l) n - m < a: n _i • 



<.***<> -^m+l ^ *£m- 



(226) 



(227) 



(228) 
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Also, 

• (-l) r_1 > 0. (229) 

If, in addition, r > 3, then 

2 ■ X'm-\-l 2 • X>f n 

(-l) r - 3 • Xm -(r- 3) > (-lY-' ■ 2 " m 3 (r - 4) , (230) 

and 

{~l) r - 2 • x m - (r - 2 ) > (-l)^ 3 • (231) 



Proo/. We obtain (|2"2"g)) from the combination of (g55P and Theorem [T31 Also, we obtain 
from the combination of Theorems I2TJ1 (see also the proof of Theorem |2"T|) . Finally, we combine 
Theorems IM1I2TJ1 and Corollary [5] to obtain (ggQD , (j2"3"Tj) (see also the proof of Theorem EH]), ■ 

In the following theorem, we estimate the accuracy to which the complex numbers otj of Theo- 
rem [T7] are evaluated. 

Theorem 26. Suppose that the n by n symmetric tridiagonal matrix A is that of Definition [7] in 
Section \4-l\ and that A is an eigenvalue of A. Suppose also that 1 < k < n and 1 < I < n — k — 1 
are integers, and that 

\-A 1 >--->X-A k >2>\- A k+1 > • • • > A — A k+l+1 > 0. (232) 

Suppose furthermore that the real numbers 9 k , • ■ • , k +i are defined via (| 130[) in Theorem \14\ Then, 

COS(^) - COS(0j+i) + i 

< 1 i ' 



sm(6 j+1 ) • sin ( flj,+ * J+1 ) - ^fc+i 



/or every integer j = k+l,...,k + l — 1. j4Zso, 

1 



l-cos(6>j) 
/or every integer j = k+l,...,k + l— 1. 
Proo/ We combine ([T3"0l of Theorem [Til with (l2"3"2"j) to obtain 



< 1 J" ' (234) 



9 J+ i) cos(0j) — cos(0j + i) ^_ cos(9j) — cos(9j + i) 



ij+i ) ■ sin 



(l-cos(6» J ))-(l + cos(6» i )) l-cos(^) 



2 — A + Aj+i ^4j+i — 
which implies (|233[) . By the same token, 

1 2 



^•+2-^+1 < ^ +2 ~^ +1 , (235) 



l-cos(^) 2-A + Aj+i /Ij+i-Afc+i 

which implies (j234[) . 



< I ^> (236) 
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Corollary 9. Under the assumptions of Theorem \2b\ suppose, in addition, that k + l<j<k + l — l 
is an integer, and that 

t +2 ~A +1 < \ (237) 
Then, the evaluation of ctj+i from ctj via (I150p in Theorem \ 1 7\ is stable. 

Proof. The statement follows from the combination of Theorem [5U and Corollary 2] ■ 

Corollary 10. Under the assumptions of Theorem \2b\ suppose, in addition, that k + 1 < j < k+l — 1 

is an integer, and that the vector (xj,Xj + i) e R 2 is given to relative accuracy e > 0. Then, the 
vector (aj,/3j) G C 2 , defined via (|147p in Theorem \ll\ has relative accuracy r(atj), where 

r{a 3 )<e-- — . (238) 

A 3 + l ~ A k+1 



Proof. The statement follows from the combination of Theorem [55] and Theorem [TU ■ 
4.3 Error Analysis for Certain Matrices 

In this section, we consider a certain subclass of tridiagonal matrices, described in Definition [1] For 
the matrices of this class, we will present a detailed analysis of some of the quantities, introduced 
in Section [OJ 

Due to Theorem [521 and Corollaries [SJ [TU1 the precision of the complex parameters ctj of The- 
orem [T7] is determined by the magnitude of the right-hand side of (|233p , (I234p in Theorem [551 
Obviously, the latter depend on the choice of the sequence < f\ < f% < . . . in Definition [T] In the 
following two theorems, we will estimate these quantities for several special cases. 

Theorem 27. Suppose that a,c,K> 1 are real numbers, that n < K 1 / 11 ■ c is an integer, and that 
Ai , . . . , An are defined via 

Aj = 2+(J^ , (239) 

for every j = l,...,n. Suppose also that the n by n symmetric tridiagonal matrix A is that of 
Definition [H that A is an eigenvalue of A, and that 

A = 4 + K ■ c- b , (240) 

for some real < b < 1. Suppose also that 

b<^, (241) 
a + 2 

that 1 < k < n and I > 3 are integers, that 

\-A 1 >--->\-A k >2>\- A k+1 > • • • > A — A k+l+2 > 0, (242) 

and that 

^± i ±i >C0S (_i_)>A^«. (243) 
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Then, 

2 



A k +i - A k 

id also 

Ak+i+i — Ak+i 



= Q 2a/(a+ 2 y 



Ak+i — Ak 



o (c-v+ b - b /")/A = o L-v^ , c 



oo. 



Proof. First, we combine (|239[) with (|240[) to conclude that 

fc" =if-C a - 6 -(l+o(l)), C^OO. 

Next, we assume that 

I = D ■ c d ■ (1 + o(l)) , c^oo, 
for some real numbers -D, d > 0, and that 

— = o(l), c — » oo. 
fc 

To validate these assumptions, we combine (|246l) . (|247l) . (|248D with (|243[) to obtain 

A - A k+ i = A - A k _ A k+i - A fc _ 
2 ~ 2 2 

We combine (|2"5§1) . (|2lrJl) and (12"4U)) to obtain 



(J)-(l + o(l)), 



c — > 00. 



(k + l) a - k a _ 7T 2 
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We combine (|246jl . (|248l) and (I250| to obtain 



r 2 „a 



ak a ~ L l = — ■ — .(l + o(l)), c^oo, 



and substitute ([MrJl) . ([M7|l into (|2"5Tj) to obtain 



4a 

We combine (|2T71) . ([2321 to obtain 



I s = — ■c a -k 1 ~ a - (l + o(l)) 
4a 

= ^-c Q -(^- C a - b )^-(l + o(l)) 

' .^+b-b/a .(1 + (1)), C ^oO. 



' 2 7-.- 1 ~° \ 3 



4a y 
and 

d _ 1 + 6-b/g 
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Due to the combination of (|241[) and (12541) . the assumptions (|247[) , (|248[) were consistent. We 
combine (|2T7|> . ([215]) . (j2"55j) . (gMP to conclude that 



A k+ i - A k it 2 

= 0(l). C 3- (1+b(1 - 1/a)) .(1 + 0(1)), c^oo. (255) 

We use (flHTj) to obtain 

2 / , ( 1\\ 2 3a 2a 

3-( 1 + 6 -( 1 -ajj < 3-^^^ (256) 

and substitute (j256)l into (l255j) to obtain ([2441 . Also, we combine ([246]) . ([247)l . (|248l) with (|239l) . 
([M31 to obtain 

A k+l+1 - A k+i = ak 11 - 1 ■ c- a • (1 + o(c)) 

= 0(1) ■ c -«+(«-&M«-l)/a . (1 + o( c )) 

= O(l) ■ c-( 1+b - b / a ) -(l + o(c)), c^oo. (257) 

Finally, we combine (g55J) , (g5g} and ([2371) to obtain (|2"4l)l) . ■ 

Theorem 28. Suppose that a, c, K > 1 and are rea? numbers, that n < K 1 / 11 ■ c is an integer, and 
that A\ , . . . , A„ are defined via 

A,-=2+Q , (258) 

/or every j — l,...,n. Suppose also that the n by n symmetric tridiagonal matrix A is that of 
Definition Ql that A is an eigenvalue of A, and that 

\ = A + K-c- b , (259) 

for some real < b < 1 . Suppose also that, as opposed to (|241l) in Theorem \21\ 

(260) 

a + 2 

that 1 < k < n and I > 3 are integers, that 

\-A 1 >--->\-A k >2>\- A k+1 >■■■> X- A k+ i +2 > 0, (261) 



and that 



Then, 



and also 



2 = O (>/(°+ 2 )) , c ^oo, (263) 



U+i — A k 



A k +i+i — A k+ i _ 



A 



k+l — ^Ife 



( c -«/(«+ 2 )) = O (c- 1 / 3 ) , c -j- oo. (264) 
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Proof. First, wc combine (|258[) with (|259[) to conclude that 

k a = K ■ c a - b ■ (1 + o(l)) , c^oo. (265) 

Next, we assume that 

I = D ■ c d ■ (1 + o(l)) , c^oo, (266) 
for some real numbers D, d > 0, and (as opposed to (|248l0 that 

j = 0(1), c^oo. (267) 

To validate these assumptions, we combine (|265l) . (|266D . (|267l) with (|262[) to obtain 

A - A k+l _ A - A k _ A k+ i - A k ^ 7T 

2 ~ 2 2 

We combine ([2351) . (|2"rJ51 and (l2"rJgf to obtain 



cos (^)- (l + o(l)), c^cx.. (268) 



(k + l) a ~k a TT 2 , , x 

" ^ = 4,2 -(1 + 0(1)). c->oo. (269) 

We combine (|265)) . (|267f and (|269| to obtain 

Z° = O(l) • ^ • (1 + o(l)) , c^oo, (270) 

and substitute (|2"rJ5l) . ([2^6]) into (j2T0")) to obtain 

? = 0(l)-c a/(Q+2) -(l + o(l)), c^oo. (271) 

Due to the combination of (|260[) and (|271l) . the assumptions (|266[) . f|267[) were consistent. We 
combine (j2"6l)j) . pjgj) . (j2"TTj) to conclude that 

2 8 • I 



_.(1 + 0(1)) = 0(c 2q /( q + 2 )), c^oc, (272) 

which implies ([2"6"5jl . Also, we combine (|2^5|) . (l2l)6"|) . (l2l)T|) with (gSBD , (|2"52"|) to obtain 
- A fc+ , = a • P- 1 • c- a • (l + o(c)) 

= O (Va+a(a-l)/(a + 2)^ = Q ^-3a/(a+2)^ ^ fi ^ ^ (373) 

Finally, we combine (l2"72l and (j2"731 to obtain (|2"M|) . ■ 

Remark 15. Suppose that a, fc > are reaZ numbers. Suppose also that the function g : (0, oo) — > K 
is defined via the formula 

/or a/Z reaZ x > 0. Then, g is monotonically decreasing. 
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Remark 16. Due to the combination of Theorem \10\ with (|258[) . (|259l) . any eigenvalue X of A 
satisfies the inequality 

A< A n + 2 = 4+(-j =4 + K = 4 + K-c°. (275) 



The following theorem is in the spirit of Theorems [571 12"51 

Theorem 29. Suppose that a,c,K> 1 are real numbers, that n < K 1 / 11 ■ c is an integer, and that 
A\ , . . . , A n are defined via 

Aj = 2+(A , (276) 

for every j = 1, ...,n. Suppose also that the n by n symmetric tridiagonal matrix A is that of 
Definition]]^ that A is an eigenvalue of A, and that 

\ = 4 + K ■ c- b , (277) 

for some real < b < 1. Suppose also that 1 < m < n and r > 3 are integers, that 

> A - A m -r > • • • > A — Am > -2 > A - A m +u (278) 



id that 



Then, 



and also 



A m A. 



>cos(^- T > . (279) 



= 0(c 2/3 V c->oo, (280) 



m -^m — r 



A — A 



(c~ 1/3 ), c^oo. (281) 



Proof. We combine (|276p and (|277D to obtain 

m = 4 1/Q • c- (1 + o(c)) , c->oo. (282) 

Next, we assume that 

r = L>-c d -(l + o(l)), c^cx), (283) 
for some real numbers D, d > 0, and that 

— = o(l), e->-oo. (284) 
m 

To validate these assumptions, we combine f|282[) . (|283p . (|284p with (|279|> to obtain 

:os(|^) -(1 + 0(1)), c^oo. (285) 
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We combine (1276) . ([252) and ([2551) to obtain 

m a — (m — r) a n 2 , , NN , 
^ L = i^ + c ^°°- ( 286 ) 

We combine (p?521) . ([2541 and (EM)) to obtain 

r ■m a - 1 = 0(1) • ^ ■ (l + o(l)), c^oo, (287) 

and substitute ([252]) . ([2531) into p57|) to obtain 

r = 0(c 1/3 ), c^oo. (288) 



Due to the combination of (|282[) and (]288[> . the assumptions (|283[) . (|284|) were consistent. We 
combine ([283"[) . (j285[) . (j288[) to conclude that 

2 — 2 



^•(l + o(l))=o(c 2 / 3 ), c^oc, (289) 

which implies (jMij) . Also, we combine (12521) . (12531) . (EMI with <P76l - (|279l to obtain 

A m - r - A m - r - X = a ■ m a ~ 1 • c~ a • (1 + o(c)) 

= of~Y c^oo. (290) 

Finally, we combine ([2551) and P§0l to obtain (|25T|) . ■ 

In the following theorem, we summarize Theorems [271 HFJ H3 

Theorem 30. Suppose that a > 1 and c > 1 are reaZ numbers, and that n > is an integer such 
that 

n = 0(c), c-^oo. (291) 
Suppose also that A\, . . . , A n are defined via 

A j= 2+(^ , (292) 

for every j — 1, ...,n. Suppose also that the n by n symmetric tridiagonal matrix A is that of 
Definition [TJ that A is an eigenvalue of A, and that 

A > 4. (293) 

If 

\-A 1 >--->\-A k >2>\- A k+1 > • • • > A — A k+l+2 > 0, (294) 



^^>co S (^)>A^±i (296) 



3G 



for some integer 1 < k < n and I > 3, then 

2 



A k+ i - A k 



O 



(c 2a /( a+2 )), c^oo, (296) 



id also 



= O (V 1/3 ) , c -> oo. (297) 



> A - A m _ r > ■ ■ • > A — A rn > — 2 > A - A m+i , (298) 



> cos ^— r ^ s ( 2 ") 



2 V 2r - 1 

/or some integer 1 < m < n and r > 3, i/ien 

2 



1 m -^-m— 7" 



(c 2/3 ) , c->-oo, (300) 



= 0( C - 1/3 ), c^oo. (301) 



A m A, 

id also 

A — A 

A — A 



Remark 17. While a more careful analysis allows to replace the asymptotic estimates (|296p . (|297p . 
pOOp . (|301[) with explicit non- asymptotic inequalities, the proofs, albeit straightforward, are somewhat 
long, and will be published at a later date. 

5 Numerical Algorithms 

In this section, we describe several numerical algorithms for the evaluation of the eigenvectors of 
certain symmetric tridiagonal matrices. 

5.1 Problem Settings 

In this subsection, we describe the problem we want to solve. 

Assumptions. Suppose that n > is an integer, that < /i < fa < . . . is monotone sequence 
of positive real numbers, and that the numbers A±,...,A n are defined via (|65p of Definition [1] 
Suppose also that the n by n symmetric tridiagonal matrix A is defined via (|64p of Definition [TJ and 
that the real number A is an eigenvalue of A. 

Goal. Our goal is to evaluate the unit-length eigenvector 

X = {X 1 ,...,X n )eW n (302) 

of A corresponding to A, whose first coordinate is positive (i.e. Xi > 0). 

Observation. We observe that A > (see Theorem O, and that A < 2 + A n (see Theorem ITO]) . 
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Equivalent formulations of the goal. Obviously, the problem of evaluating x given A and A 
is equivalent to computing the solution of the linear system 



B ■ X = 0, (303) 

where B = A—X-I. Also, due to Theorem[7J this problem is equivalent to evaluating the solution to a 
certain three-terms recurrence relation (namely, the one specified via (|66[) . (|67[) . (|68[) of Theorem [7]). 

Desired numerical properties of the solution. We want our solution to possess the following 
property: each coordinate Xj of the vector X should be evaluated with high relative precision. In 
other words, suppose that X = (X±, . . . , X n ) £ 1" is our numerical approximation to X. For every 
j = 1, . . . , n, we want the absolute value of the quantity Ej, defined via 

e 3 = (304) 

to be small. 

Observation. This task is potentially difficult if \Xj\ is small compared to ||X|| = 1. For exam- 
ple, if |Xl| < e, where e is the machine precision (e.g. e ~ 10~ 16 for double-precision calculations), 
it is not obvious why |ei| should be less than 1, let alone be "small" (see also Section [JJ). 

Relation between Xj-t, Xj and Xj+i- For every j = 2, . . . , n — 1, the relation between three 
consecutive coordinates Xj-i, Xj and Xj + \ of the vector X is expressed via (|67[) of Theorem [71 
more specifically, 

AVi + {Aj - A) • Xj + X j+1 = 0, (305) 

for every j = 2, . . . , n — 1. It turns out that the qualitative behavior of Xj—i, Xj and Xj + \ relative 
to each other depends on A — Aj in the following way. If A — Aj > 2, then all the three coordinates 
have the same sign, and < \Xj\ < \Xj+i\ (see Theorem [TTj) . If X — Aj < —2, then the signs of 

Xj-i,Xj+i are opposite to the sign of Xj, and > \Xj\ > \Xj+i\ (see Theorem[JJ]). Finally, if 

—2 < X — Aj < 2, then the relation is somewhat more complicated (see, for example, Theorems [T71 

EH Oil 125). 

Assumption on A. In the view of the latter observation, we will consider the case in which the 
coordinates of X exhibit all the behaviors described above (that is, in this sense, the most general 
case). This is achieved by making the following assumption on A. Suppose that 

l<k<k + l + l<p<p+l<m-r<m<n (306) 

are integers. Suppose also that 

A-^i>--->A- J 4 fe > 2> 
A - A k+1 > • • • > A - A k+l+1 >--->X-A p > 0> 
A - A p+1 >■■■> X - A m - r > ■ ■ ■ > X — A m > -2 > 

A - A m+ i > ■ ■ ■ > X — An. (307) 



Suppose furthermore that 



and that 



A - A k+l+ i ( 7T \ A - A k+ i +2 . . 
2 >C ° S (21^1 ~ 2 ' (308) 



A m -r — X f 7T \ A m _ r _i — A fonc\\ 

^^ >COS 273lP 2 • ( 309 ) 
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Observation. We combine p07[) with Definition [JJ to conclude that 



4 + /i<A</„. (310) 

If the left-hand side inequality in (|310[) did not hold, then the region of growth would be empty. If 
the right-hand side inequality in piOl) did not hold, then the region of decay would be empty. While 
the obvious simplification of the algorithm of the next section will handle such cases, we discuss only 
the general case for the sake of simplicity of presentation. 

5.2 Informal Description of the Algorithm 

This section contains an informal description of an algorithm for the evaluation of X = (X\ , . . . , X n ) G 
W 1 (see (|302p V In Section [5731 the algorithm is described in a more precise and detailed way. Sec- 
tion [ITU contains an outline of the steps of the algorithm. 

For any A— eigenvector x = (x\ , . . . , x n ) € K™ of A and every j = 2, ...,n— 1, the three 
consecutive coordinates Xj-i, Xj, Xj + i satisfy the identity (|S7|) of Theorem [7] (see also (|305[) above). 
If A — Aj > 2, we say that j is in the region of growth. If A — Aj < —2, we say that j is in the region 
of decay. Otherwise, we say that j is in the oscillatory region (see Definition [2]) . 

First, we evaluate the coordinates of x in the region of growth by setting the first coordinate to 
be equal to 1 and "going forward" (i.e., Xj+i is computed from Xj and Xj^i via (1571) ). We use the 
same recurrence relation to evaluate the first few coordinate in the leftmost part of the oscillatory 
region (when A — Aj is still close to 2). All of these coordinates are positive, and this computation 
is numerically stable (see Sections 15.3.11 15.3.21 below for more details) . 

Then, up to the middle of the oscillatory region (i.e. up to j = p with A — A p w 0), we 
"complexify" the coordinates of x by relating to each pair (xj,Xj+\) a complex number aj, whose 
real part is Xj/2. These a/s are evaluated as follows. The "leftmost" aj is evaluated from the 
corresponding Xj,Xj+\ via solving a certain two by two linear system; in this computation we lose 
several digits (see, for example, Corollary [T0| . However, the evaluation of each subsequent aj from 
its predecessor is a stable procedure (see, for example, Corollary[S]and Section l5.3.3l for more details). 
Once we have all ay's, each Xj in the left half of the oscillatory region is evaluated as 2 • $t(aj) (see 
Section 15.3.31 for more details). 

At this point, we have evaluated the left coordinates x±, . . . , x p , Xp+i of the A— eigenvector x of 
A, whose first coordinate is x\ = 1. 

Suppose now that x — {x\, . . . , £„) € R™ is the A— eigenvector £ of A whose last coordinate is 
x n = 1. The evaluation of the right coordinates of x mirrors the evaluation of the left coordinates of 
x, with the main difference being that the direction is reversed. More specifically, we set x n = 1, and 
evaluate the coordinates of x in the region of decay via the recurrence relation (p7|) in the backwards 
direction (i.e. ij-\ is evaluated from ij and ij+i). We use the same recurrence relation to evaluate 
the last few coordinate in the rightmost part of the oscillatory region (when A — Aj is still close to 
—2). All of these coordinates have alternating signs, and this computation is numerically stable (see 
Sections 15.3.41 15.3.51 below for more details) . 

The evaluation of the rest of the coordinates of x in the right part of the oscillatory region (i.e. 
for > A — Aj ) is similar to the evaluation of the coordinates of x in the left part of the oscillatory 
region (see also Section 15.3.61 below for more details) . By the end of this computation, we have the 
coordinates i p , x p +\, . . . , x n of i. 

Since the A— eigenspace of A is one-dimensional (see Theorem[7J, x and £ must be, up to a scaling, 
the coordinates of the same eigenvector (in exact arithmetics). We use this observation to rescale 
x p , . . . ,x n in such a way that (x p , x p +i) = (x p , x p +\). By doing so, we obtain all the coordinates of 
the A— eigenvector x of A whose first coordinate is x\ ~ 1. We normalize x (i.e. divide it by ||x||) 
to obtain the unit-length eigenvector X whose first coordinate is positive (see Sections 15.3.71 15.3.81 
below for more details). 
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5.3 Detailed Description of the Algorithm 

We evaluate X = (X 1 ,..., X n ) € R" (see (gD5)> ) as follows. 



5.3.1 The region of growth: x\,..., Xk+i 

First, we set x\ = 1, and evaluate X2, ■ ■ ■ , %k+i iteratively via (JHSJ), ((57)) of Theorem [7] These 
coordinates are all positive, and grow as the indices increase (see discussion in Section [5. 11 and also 
Theorem 1111 for more details). Also, this computation is stable numerically (see Theorem 1121 for 
more details). Intuitively, it says that three-term recurrences are stable in the growing direction. 

5.3.2 The leftmost part of the oscillatory region: Xk+2, %k+i-i 

We observe that < A — Ak+i < 2, due to (|307[) . In other words, we have reached the oscillatory 
region (see Definition [2]) . In particular, the coordinates are not expected to have the same sign, 
the monotonicity is not guaranteed, and, moreover, some coordinates might be significantly smaller 
than the other (for an extreme example, see Theorem 2]). In other words, the use of the recurrence 
relation (p7|) to evaluate the coordinates in the oscillatory region can potentially result in loss of 
accuracy. 

On the other hand, the first few coordinates in the oscillatory regime do not change too fast, 
and are all positive (see Theorem 1211 for a more precise statement). We compute the coordinates 
Xfc + 2, ■ • ■ j^fc+i-i iteratively via (|67l) . moreover, this computation is numerically stable (see Theo- 
rem [22]). We note that I is determined from (|308l) : in particular, A — Ak+i+i is still fairly close to 
2. 

5.3.3 Up to the middle of the oscillatory region: Xk+i, ■ ■ ■ , x p +\ 

Up to the middle of the oscillatory region (A — A p w 0), the coordinates of x can vary in both sign 
and order of magnitude (see, for example, the eigenvectors of the matrix in Theorem [4]) • Therefore, 
if we still use ([67]) in this region to evaluate the coordinates one by one, "small" coordinates might 
be computed as a sum of two "larger" coordinates, resulting in a loss of accuracy; moreover, the 
stability properties of such computation are not obvious. 

To bypass these obstacles, we complexify the coordinates in the following way. First, we com- 
pute the complex number ak+1-2 from Xk+1-2, Xk+1-1 via the two by two linear system (|147[) of 
Theorem [T7] (the solution is obtained via Q137I) of Theorem fT5T) . The condition matrix of the corre- 
sponding two by two linear system is roughly I 2 , where I is defined via (|308[) (see Theorem[T(3 (|130[) 
of Theorem 1 141 and ([308])). In other words, in the evaluation of ak+1-2 from Xk+1-2 and x^+i-x we 
expect to lose 2 -log 10 (Z) decimal digits (see e.g. Corollary [TU] also, for an example for this condition 
number for a certain class of matrices A, see (|296[) of Theorem [30]) . 

We then evaluate a^+i-i, ■ • ■ , ot P -i, a p iteratively via (|150[) of Theorem [T71 In other words, 
each oij+i is obtained from ay via rotation by the angle 9j and adding a pure imaginary correction 
term. This computation is stable provided that the correction term is small compared to |oy| (see 
Corollary [H Theorem 121)1 and Corollary ^ . In Section |4~51 we demonstrate that, for a certain class 
of matrices, the correction term is indeed small, and thus all ak+i-2, ■ ■ ■ ,a p are evaluated with the 
same relative accuracy (see Theorem [53] Remark fTS] and Theorem I3TJ1) . 

Next, we evaluate Xk+i, ■ ■ ■ , x p -i,x p via fll48[) of Theorem [T7l In other words, Xj is the real part 
of aj/2, for each j = k + l, . . . ,p. In particular, each Xj is computed with the same absolute accuracy 
as ay. On the other hand, each Xj is expected be evaluated to roughly log 10 (|aj/a;j|) correct decimal 
digits less than ay. 

Finally, we evaluate x p +\ from a p via (|149[) of Theorem [T71 
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5.3.4 The region of decay: x m , ...,x n 

The rest of the coordinates of x are computed in the reverse direction, i.e. we iterate backwards from 
n to p. We will denote these coordinates by x p , . . . , x n , to distinguish them from X\, . . . , 
whose evaluation is described in Sections 15.3.11 15.3.21 15.3.31 

First, we set x n = 1, and evaluate x n -\, . . . ,x m iteratively, in the backward direction, via (|68p . 
(|57|) of Theorem [71 All x m , . . . ,x n have alternating signs, and their absolute values decay as the 
indices increase (see and Theorem [13] for more details). Also, this computation is stable numerically 
(see Theorems [13J [12] for more details) . Intuitively, it says that three-term recurrences are stable in 
the growing direction. 

We note that the evaluation of x m , . . . , x n mirrors the evaluation of x\, . . . , Xk+i, described in 
Section [5.3.11 above . We also observe that, depending on the parity of n — m, the coordinate x m can 
be either positive or negative. If x m < (that is, n — m is odd), we flip the signs of all x m , . . . , x n . 
Obviously, the recurrence relation ([571) . (|55| is still satisfied by x m , . . . ,x n ; on the other hand, we 
have made sure that x rn is positive. 

5.3.5 The rightmost part of the oscillatory region: x m _/ r -2)t ■ ■ ■ , %m-x 

As we enter the oscillatory region from the right, the magnitudes of the first few x/s do not change 
too fast, and they all have alternating signs (see Theorem [231 for a more precise statement). We 
compute the coordinates x m _t r -2)> ■ ■ ■ i %m-i iteratively via (1571) of Theorem [7J in the backward 
direction. Moreover, this computation is numerically stable (see Theorem [2"2"jl . We note that r is 
determined from (I309|) : in particular, A m - r — A is still fairly close to 2. 

We note that the evaluation of £ m _( r _2)> • • • i im-i mirrors the evaluation of Xk+2, ■ ■ ■ , Xk+l—li 
described in Section 15.3.21 above. 

5.3.6 Down to the middle of the oscillatory region: x p , . . . , i m ( r i) 

The evaluation of x p , . . . , x m _(, — i) mirrors the evaluation of Xk+i+i, ■ ■ • , £p+i, described in Sec- 
tion [533] above. In other words, we complexify ij's by introducing 7j's (the real part of each 7 3 is 
twice Xj), and evaluate 7j's iteratively, in the backward direction. 

More specifically, we first compute 7 m -( r _i) from £ m _( r -2)j i m -(r-3) by solving the two by two 
linear system (j218p of Theorem [24] (the solution is obtained via (11 3T[) of Theorem [TBI . The condition 
matrix of the corresponding two by two linear system is roughly r 2 , where r is defined via (|309p (see 
Theorem [T51 (j213[) of Theorem [231 and (J30S]))- in other words, in the evaluation of 7 m -( r -i) from 
Xm-(r~2) an d x m _( r _ 3 ) we expect to lose 2 • log 10 (r) decimal digits (see e.g. Corollary [TU] also, for 
an example for this condition number for a certain class of matrices A, see (|296p of Theorem I30|) . 

Next, we evaluate 7 m _ r , . . . ,7 P ,7 P _i iteratively, in the backwards direction, via (j22ip of The- 
orem [2"4l This computation is similar to the evaluation of otk+i-i, ■ ■ ■ , a^-i, ct p , described in Sec- 
tion 15.3.31 and has similar properties in terms of accuracy and stability (see also Theorem 1261 
Corollary [9] Remark [TBI and Theorem l30l) . 

Then, We evaluate %rn — (r—l) i %m~ri • ■ • ? 

from 7 m _( r+ i), . . . ,7 P -i, respectively, via (j219p of 

Theorem 1231 

Finally, we evaluate x p from 7 P -i via (I220p of Theorem [231 

5.3.7 Gluing xi, . . . ,x p , x p+ i and x p , x p+ i, . . . , £„ together 

We have evaluated the first p + 1 coordinates of the A— eigenvector x of A whose first coordinate is 
xi = 1. Also, we have evaluated the last n — p+1 coordinates of the A— eigenvector x of A whose last 
coordinate is x n = (— l)" _m . However, due to Theorem[7J the A— eigenspace of A is one-dimensional. 
We use this observation to conclude that, up to scaling, x p , . . . , x n would be, in exact arithmetics, 
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the last n — p + 1 coordinates of the eigenvector x. In other words, in exact arithmetics, if the real 
number s is defined via 

^ = {x p /x p , if \x p \ > \x p+1 \ ^ 0, 

1 Xp+i/xp+i, otherwise , 

then the vector 

[x\ , . . . , X p: X p +i , S • Xp-j-2 j • • • : S ' ^n) (312) 

is the A— eigenvector of A whose first coordinate is x\ = 1. We observe that s ^ 0, since Xp + xp+i 
is strictly positive, due to Theorem [7] Thus, we evaluate x p +2, ■ . ■ ,x n via the formula 

Xj — 5 ■ x j j ^313^ 

for every j = p + 2, . . . , n, where s is defined via pill) . 

We observe that s is evaluated from x p , x p +i, x p , x p +i, and each individual in the oscillatory 
region might be evaluated less accurately than aj (see Section l5.3.3[) . Nevertheless, we observe that 
9 p w 7r/2 and <p p -i ~ tt/2, due to the combination of Theorems [TU [531 with (I307|) . We combine 
this observation with ([Tig]) . ([13511 of Theorem [T71 (|2T^1) . (|2"2"0|l of Theorem [Ml Corollary [3] and 
Corollary \8\ to conclude that 



l^l + K+il 



4K| 2 , 



|2 + |x p+1 | 2 «4| 7 p^ 1 |. (314) 



It follows from the combination of (I311[) and ([31411 that s is evalauted to roughly the same relative 
accuracy as a p /"f p -i. 

5.3.8 Normalization 

We define the real number d via the formula 



d= x /xj + --- + x 2 n . (315) 

In exact arithmetics, the vector 



X 



(X 1 ,...,X n )=(^,...,^) (316) 



is the unit-length A— eigenvector of A, whose first coordinate is positive. 

In Section 15.3.31 we mentioned that, in the oscillatory region, the coordinates Xj might have a 
significantly lower relative accuracy than the corresponding ay. Motivated by this observation, we 
evaluate d without using the coordinates Xj in the oscillatory region explicitly, in the following way. 

First, we evaluate the real numbers di via the formula 



v 



4 (l^| 2 +cos(^-)-5R(a^)) 
j=k+l-i 
p 

4 |a J f-(l + cos(6' J )-cos(6l J +2Arg(a J ))), (317) 

j=k+l-i 
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where Arg(z) = £ for any complex number z = R ■ e*^ . Due to Corollary [51 



p 



= |x fe+; _ 1 | 2 + |x p+1 | 2 + 2 \ x 3?- ( 318 ) 



Also, we evaluate the real number d r via the formula 



Due to Corollary [51 



m— (t— 1) 

4 E (l7i| 2 +co S fe)-5i(7^)) 
i=p 

(r — 1) 

4 E l7 J | 2 -(l + cos(^-)-cos(^+2Arg( 7j ))). (319) 



m-(r— 2) 

d r = |^ +1 | 2 + |x m _ (r _ 3) | 2 + 2 E l^l 2 - ( 32 °) 



We combine (|3TT|) . (pJTgj) . (|52D| to obtain 



m — (r — 3) 

d ( +s-d r = |a; fc+; _ 1 | 2 + |a; m _ (r _3 ) | 2 + +2 ^ W 2 - ( 321 ) 



Therefore, in exact arithmetics, 

d 2 = 

I 12 I I |2 I j I j /c+Z — 2 n 

+ l g "-fr-3)l + J- \ Xj \ 2 + J2 N"' ( 322 ) 

j — 1 j= m— r+4 

We compute d by evaluating the right-hand side of (I322p numerically. 

The relative accuracy to which d is evaluated will be dominated by the relative accuracy of the 
summands in (|317l) . (|319l) (each one of these summands is obviously positive). 

Finally, we obtain the unit-length eigenvector X by normalizing the eigenvector x. This is done 
by dividing every coordinate of x by d (see (|316p ). 



5.4 Short Description of the Algorithm 
Step A: evaluation of the left coordinates. 

1. Set xi = 1. 

2. Compute X2 via ([66)) of Theorem [3 

3. Compute £3, ... , x^, Xk+i iteratively via (pTj) of Theorem [7J 

4. Compute Xk+2, ■ ■ ■ , £fc+;-2, iteratively via ([o7|) of Theorcm[7J 

5. Compute afe+/_2 from Xk+1-2, Xk+1-1 via (|147l) of Theorem [171 

6. Compute ak+1-1, ■ ■ ■ , a p _i, a p iteratively via (|150l) of Theorem IT71 

7. Compute a^+i, . . . , x p -\,x p via (|148l) of Theorem fTTl 

8. Compute a; p+ i from a p via (|149l) of Theorem [171 



Step B: evaluation of the right coordinates. 



43 



1. Set x n = 1. 

2. Compute £„_i via of Theorem!?] 

3. Compute x n - 2 , • ■ • , x m +i, x m iteratively via (1571) of Theorem [7J 
If i m < 0, flip the signs of all x m , . . . , x n . 

4. Compute x m -i, ■ ■ ■ ,^ m -(r-2) iteratively via (J57J of Theorem)?) 

5. Compute j m _r r _i\ from i m _( r -2),im-(r-3) v i a Q218[l of Theorem [Ml 

6. Compute 7 m _ r , . . . , 7 P , 7 P -i iteratively via ()221[) of Theorem [Ml 

7. Compute f m _( r _i-j, x m - r , . . . , x p +i via (|219[) of Theorem [Ml 

8. Compute x p from 7 P _i via (j220[) of Theorem [Ml 



Step C: glue them together. 

1. If x p • x p + x p+ i ■ x p+ i < 0, flip the signs of all x p , . . . , x n . 

2. Set s —— Xpjxp. 

3. If |x p +i| > |x p |, set s = Xp+i/ip+i. 

4. For every j = p + 2, . . . , n, compute Xj via the formula 

5. Compute d = \J x\ + ■ ■ ■ + x\, via evaluating the right-hand side of (I322p . 

6. For every j = 1, . . . , n, compute Xj via Xj = Xj/d. 



5.5 Error Analysis 

In Sections 15.31 15.41 we described an algorithm for the evaluation of the unit length A— eigenvector 
X = (Xi, . . . ,X n ) of A, whose first coordinate is positive. In this section, we analyze the relative 
accuracy to which each coordinate is evaluated. In this analysis, we use the following notation. 
Suppose that y is a quantity to be evaluated, and eval(y) is a numerical approximation to this 
quantity. We define the absolute error abs(y) via the formula 

abs(y) = \y- eval(y)\. (323) 

Also, we define the relative error rel{y) via the formula 

y - eval(y) 



rel(y) = 



y 



(324) 



Loosely speaking, — log 10 (reZ(y)) is the number of decimal digits in which y and eval(y) coincide. 
Our analysis is based on the following observations. 

Observation 1. All of x\, . . . , Xk, Xk+i have roughly the same relative accuracy (see Sec- 
tion [53H1)- This is despite the fact that X\ can be significantly smaller than Xk (see Theorem ITTj) . 

Observation 2. Also, all of Xk, £k+li • • • i x k+i-i have roughly the same relative accuracy (see 
Section EH}. 

Observation 3. All of x m ,x m+ i, . . . ,x n have roughly the same relative accuracy (see Sec- 
tion [5331) ■ 

Observation 4. Also, all of x m - r +2, ■ ■ ■ , x m , x m +i have roughly the same relative accuracy (see 
Section 15X5)1 . 

Observation 5. The relative accuracy to which ak+1-2 is evaluated is roughly 

rel(a k+ i- 2 ) ~ rel(x k+ i^ 2 ) • -j — j- ~ rel(x k+l - 2 ) • I 2 (325) 

A k +l — Ak 

(see Section 15. 3.31 and also Theorems 1T51 [2l)l Corollary ITU)) . 
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Observation 6. The relative accuracy to which ak+i-i, ctk+i, • ■ • , a p are evaluated is roughly 
the same, provided that, for example, 



Aj+2 - Aj+i , 1 



< (326) 



for every j = k + I — 1 , . . . , p (see Section 15.3.31 and also Theorems [T7J [Ml an d Corollaries [H ^ . In 
other words, 

rel(ctj) & rel(ctk+i-2), (327) 



for every j = k + I — 1, . . . ,p, provided that (|326[) holds. 

Observation 7. The relative accuracy to which 7 m _( r _ 1 ) is evaluated is roughly 

r-e/(7 TO _ (r _ 1) ) fa re/(x m _ (r _ 3 )) • - pa reZ(x m _ (r _ 3) ) • r 2 (328) 

^m— (r— 1) 

(see Section 15.3.61 and also Theorems [TJoJ [Ml Corollary [TU)) . 

Observation 8. The relative accuracy to which 7 P -i, 7 P , . . . , 7m-r are evaluated is roughly the 
same, provided that, for example, 

A i +2 ~f i+1 < l (329) 

for every j = p— 1, . . . , m — (r — 1) (see Section [5.3.61 and also Theorems l24l l26l and Corollaries [H 
[9]). In other words, 

rel(jj) ps rel^m-r), (330) 

for every j — p — 1, . . . , m — r, provided that p29p holds. 

Observation 9. The scaling parameter s, defined via (|31ip in Section l"5. 3.71 is evaluated with 
relative accuracy 

rel(s) ps rel(a v ) + rel{^ p ^i) (331) 

(see, for example, (|311[) . (|314|) in Section 15. 3 .7[) . 

Observation 10. The relative error of each summand in pi7p is roughly bounded by 

rel (K | 2 • (008(0,) • cos (2 • Arg^-) + 6,) + 1)) < (332) 

for every j = k + 1 — 2, . . . , p. We observe that, for any positive numbers p\ , . . . , p n and real numbers 
£!,...,£„ e (-1,1), 



Pl ■ (1 + £i) j \-p n - (1 + £») _ 

Pi H hPn 



< max{|ei|, . . . , |e„|} , (333) 



and combine (|333[) . pi7p . p25p . Q327P and p32p to conclude that the relative error of di is roughly 
bounded by 

rel(di) < rel{a k+ i- 2 ) ■ -. j- ~ rel(x k+ i^ 2 ) • I 4 - (334) 

Ak+l-2 — 
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Observation 11. By the same token as in Observation 10, we conclude that the relative error 
of d r , defined via (|319[) , is roughly bounded by 



rel(d r ) < ( ' 7 , m ( ' — , » rel(x m _ {r _ 2) ) ■ r 4 . (335) 

1 - COS(^ m _( r _i)) 



Observation 12. Suppose that e is the relative accuracy to which the quantities A — . . . , A — 
A n are given. In other words, suppose that 

rel{\ - Aj) w e, (336) 

for every j = 1, . . . ,n. (For example, if A — Aj are given to full precision, e ~ 10~ 16 for double 
precision calculations and e pa 10~ 35 for extended precision calculations.) We combine (|322p . Ob- 
servations 1-4 above, (I334[) . (|335l) and (|336[) to conclude that the real number d, evaluated via the 
right-hand side of (|322[) . has relative accuracy roughly 

rel(d) < e ■ (l 4 + r 4 ) . (337) 

We make the following conclusions from Observations 1-12. 

Conclusion 1. Suppose that the real number e > is that of (|336p . and the integers k, l,p, r, m, n 
are those of ([5D7|) . (|505j) . ([509")l . We combine ([3TB)) . ([322)1 . Observations 1-4 and ([537)1 to conclude 
that 

reZ(Xj) < e- (Z 4 + r 4 ) , (338) 

for every j = l,...,fc + / — 1 and every j = m — r + 2, . . . , n. In other word, the coordinates in the 
region of growth and the coordinates in the region of decay are evaluated to the relative accuracy, 
which is larger than the machine accuracy by only a factor of at most I 4 + r 4 . 

Conclusion 2. We combine ([3"25)l . ([527)1 . ([537)1 . ([Tig)) of Theorem Q7J and (|5T5)l to obtain 

a&s(X,-) < £ • | ay | ■ (re/(aj) + rel(d)) < e ■ |ay| • (; 4 + r 4 ) , (339) 

for every j = k + I, . . . ,p. We combine (|148[) with p39p to obtain 

rd ^^ c -wh' { '' +r ' h <340) 

for every j = k + I, . . . ,p. 

Conclusion 3. Also, We combine (|328)l . ([330)) . (|557)l . (I2T9)) of Theorem l24l and (I3T6)) to obtain 

dbsiXj) < e ■ | 7i _ 2 | • (re/( 7i _ 2 ) + reZ(d)) < E • | 7j -_ a | • (l 4 + r 4 ) , (341) 

for every j — p + 1, . . . ,r — m + 1. We combine (|219D with p4ip to obtain 

rem - £ 'wA\' {l4+r4h (342) 

for every j = j>+l,...,r — m + 1. 

Example. Suppose that a, c > 1 are real numbers, and that the n by n matrix A is that of 
Theorem [301 In other words, the diagonal entries of A are defined via (|292|) . Suppose also that 
A > 4, and that ([507)) holds. We combine (|555)l with Theorem I3U1 to obtain 

rel(Xj) <eO (V a /( a + 2 )) , c ^ oo. (343) 

for every j — 1, . . . , k + I — 1 and every j = m — r + 2, . . . , n. In other words, if, for example, 
a = 2, then the first several coordinates X±, X2, ... of the eigenvector are expected to be evaluated 
correctly to all but the last 2 • log 10 (c) digits (see also Theorem [1] in Section [1}. 
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Remark 18. Extensive numerical experiments seem to indicate that the estimates (|338[) . (|339[) . 
(|340p . (|34ip . (|342p are somewhat pessimistic. In other words, in practice the relative error tends to 
be smaller than our estimates suggest (see also Section^). 

Remark 19. It is somewhat surprising that, according to l|338p . the relative error of, say, X\ seems 
to be independent of the order of magnitude of X\. In particular, while X\ can be fairly small (see 
e.g. Theorem \ll\ and Corollary^, it still will be evaluated to reasonable relative precision. 

5.6 Related Algorithms 

In Section 1573"! 15 .4[ we presented an algorithm to solve the problem, described in Section I5TT1 In this 
section, we describe several related algorithms to solve the same problem, and briefly discuss their 
performance. 

5.6.1 A Simplified Algorithm 

The following algorithm is based on Theorem [71 and avoids complexification of the coordinates of 
the eigenvector (see Sections 15.3.31 l5T3~6|) . 

First, we set a; i = 1, and evaluate X2, x p , x p +\ iterativelv via (|66p . (|67p of Theorem!?! In other 
words, we compute all the coordinates in the region of growth and up to the middle of the oscillatory 
region via the three-term recurrence relation (1671) , in the forward direction (see Section 15.3. ip . 

Then, we set x n = 1, and evaluate x n -\, x n ~2, ■ ■ ■ , £p+ii %p via (|68p . (1671) of Theorem[71 in the 
backward direction. In other words, we compute all the coordinates in the region of decay and down 
to the middle of the oscillatory region via the three-term recurrence relation (1671) , in the backward 
direction (see Section l"5.3.4p . 

Next, we glue x\, . . . , x p , x p +i with x pi . . . , x n via the procedure, described in Section [5.3.7l 
In other words, we defined the real number s via (I31ip . and defined Xj via (|313p . for every j = 
p + 2, . . . , n. As in Section 15.3.71 we observe that (at least in exact arithmetics), the vector 

x = (xi, . . . ,x n ) £ R" (344) 

is the A— eigenvector of A, whose first coordinate is x\ = 1. 

Finally, we evaluate the unit-length eigenvector X — {X\, . . . , X n ) via the formula 

. .,X n ) = 1 ■ (x u . . .,!„) . (345) 

In other words, the normalization factor is computed by summing Xj's directly, as opposed to using 
(12221) of Section [OH 

Remark 20. Up to our best knowledge, the error analysis (similar to that of Section \ 5.5\l for 
the simplified algorithm, described in this section, is not available in the literature. Nevertheless, 
extensive numerical experiments seem to indicate that, in practice, the accuracy of the simplified 
algorithm is similar to that described in Sections \5.3l \5.4\ (see also Section^. 

Remark 21. When the coordinates of the eigenvector are evaluated via the three-terms recurrence 
(|67p . the choice of direction plays a crucial role. Roughly speaking, this recurrence is unstable in the 
backward direction in the region of growth, and is unstable in the forward direction in the region of 
decay (see also Sections ] 5. 3. 11 \5.3.4\ l- The use of this recurrence relation in the "wrong" direction 
leads to a disastrous loss of accuracy. 
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5.6.2 Inverse Power 



The unit-length A— eigenvector X of A can be obtained via Inverse Power Method (see Section 13.4. II 
for more details). This method is iterative, and, on each iteration, the approximation x^ k+1 ^ of X 
is obtained from a;W via solving the linear system 

(X ■ I - A) ■ x {k+1) = x {k \ (346) 

and normalizing the solution (i.e. we divide it by its norm). We observe that this method also 
evaluates A, but we have assumed that A has already been evaluated, by Inverse Power Method, 
Sturm Bisection (see Section F3.4.2I) . or some other means. On each iteration, we solve the linear 
system (|346p by Gaussian elimination (we recall that A is tridiagonal, and thus one iteration costs 
0(n) operations; see Remark [SJ. 

The following conjecture about the accuracy of Inverse Power Method is substantiated by exten- 
sive numerical experiments (see Section [7]) . 

Conjecture 1. Suppose that e > is the machine precision (e.g. e ~ 1CP 16 for double-precision 
calculations), and that the eigenvalue X of A is given up to e, i.e. 

rel(X) « e, (347) 

where rel is the relative error, defined via (|324p in Section 15.51 Suppose also that X — A\ > 2 ( see 
(|307[) ). Suppose furthermore that K > is an integer, and that 

K > + 1, (348) 

log(e) 

where X = {X\, . . . ,X n ) € K™ is the unit-length X— eigenvector of A. Then, after K iterations of 
Inverse Power Method, X\ is evaluated with high relative accuracy. More specifically, this relative 
accuracy is of the same order of magnitude as for the algorithms, described in Sections \5.2l\5.3\. \5.4\ 

\EM 

Remark 22. The statement of Conjecture]^ also holds for all the coordinates of X in the region of 
growth and the region of decay. In other words, suppose that k,m are those of (|306p . (|307p . Then, 
each of X\, . . . , X^+i and X m , . . . ,X n is evaluated with high relative accuracy by K iterations of 
Inverse Power Method. 

Remark 23. The ineguality (|348[) reflects on the fact that each iteration of Inverse Power Method 
can reduce the coordinates of the approximation x^ by the factor at most e _1 . In other words, if 
X\ « 1CP 50 , and, in the initial approximation, x^ = 0(1), then x^ will already be of the same 
order of magnitude as X\, and x^ will approximate Xi to a high relative precision. 



5.6.3 Jacobi Rotations 

In the view of Sections 15.6.11 15.6.21 one might suspect that virtually any standard algorithm would 
accurately solve the problem, introduced in Section [5.1l In other words, one might suspect that the 
small coordinates of X in the region of growth and the region of decay will be evaluated with high 
relative accuracy by any reasonable algorithm that computes eigenvectors. 

Unfortunately, this is emphatically not the case, and the accuracy of the result strongly depends 
on the choice of the algorithm. Consider, for example, the popular Jacobi Rotations algorithm 
for the evaluation of the eigenvalues and eigenvectors of A (see, for example, [3], [S], [S], |20j). 
This algorithm is known for its simplicity and stability, and, indeed, it typically evaluates all the 
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eigenvalues of A fairly accurately. Moreover, the corresponding eigenvectors are evaluated to high 
relative accuracy, in the sense that 



\\X - eval{X)\\ 
11*11 



\\X -eval(X)\\ w e, 



(349) 



where X is the unit-length eigenvector (see (|302|l ). eval(X) is its numerical approximation produced 
by Jacobi Rotations, and e is the machine precision. However, the coordinates of X are typically 
evaluated only to high absolute accuracy, e.g. 



where abs is defined via (|323[) . On the other hands, the relative accuracy of small coordinates will 
typically be poor. In particular, if, for example, X\ ~ 10 -50 , its numerical approximation, produced 
by Jacobi Rotations, will usually have no correct digits at all! 

These observations about Jacobi Rotations applied to the problem of Section 15.11 are supported 
by extensive numerical evidence. 

5.6.4 Gaussian Elimination 

Another possible method to evaluate X would be to solve the linear system 



by means of Gaussian Elimination (see, for example, [5], [B], [TO], [20] )■ Unfortunately, this method, 
in general, fails to evaluate the small coordinates of X with high relative accuracy (see, however, 
Section [5 .6. 2 [ where Gaussian Elimination is used several times, as a step of Inverse Power Method). 

6 Applications 

In this section, we will describe some applications of the algorithm, introduced in Section [SJ to other 
computational problems. 

6.1 Bessel Functions 

The Bessel functions of the first kind Jq,J±i, J±2, ■ • • , are defined via (|28p. (1291 in Section [XU A 
numerical algorithm for the evaluation of Jq(x), . . . , J n {x) for a given real number x > and a given 
integer n > is described in Section l3.4.3l In particular, according to Remark[S]in Section l3.4.31 the 
values Jo(x), . . . J n {x) can be obtained as coordinates of the unit length A— eigenvector of a certain 
symmetric tridiagonal matrix A(x) (see ((tJT|) . (p2"j) . (1531 ). 

The matrix A(x) belongs to the class of matrices, introduced in Definition [T] More specifically, 
the diagonal entries of A(x) are those of the matrix A of Theorem l30l with a = 1 and c — x/2 (see 



In other words, the principal algorithm of this paper (see Sections 15.21 15.31 15. 4[) can be used 
to evaluate the Bessel functions Jq, . , . , J n at a given point. Even more so, the accuracy of this 
evaluation is analyzed in Theorem [501 and Section 15751 Finally, in this case, the simplified algorithm 
(see Section I5.6.ip coincides with the well-known algorithm, described in Section 13.4.31 

See Section [7] for the related numerical experiments. 



abs(Xi) s, 



(350) 



(X-I-A)-X = 0, 



(351) 
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6.2 Prolate Spheroidal Wave Functions 



For any real number c > 0, the prolate spheroidal wave functions (PSWFs) of band limit c are 
defined in Section 13.31 The numerical algorithm for the evaluation of PSWFs, briefly described in 
Section 13 . 31 is based on computing the unit-length eigenvectors of certain symmetric tridiagonal 
matrices (namely, truncated versions of the matrices A c ' even , A c ' odd , defined, respectively, via (|41l) . 
dHJ) in Section [331). 

Strictly speaking, these matrices do not belong to the class of matrices, introduced in Definition]!] 
since their non-zero off-diagonal entries are not equal to one (see (HO])). Nevertheless, in the notation 

of gOD, 



A (c) _ Ac) 



1 + 



(352) 



for every k = 0, 1, 2, . . . , and 



<i = 7f 




(353) 



for every k = 0, 1, 2, In other words, for example, the matrix A c ' odd can be viewed as a small 

perturbation of the symmetric tridiagonal matrix A, defined via the formula 



Mo i \ 

1 A x 1 

1 A 2 1 

V '•■ '■) 



(354) 



where j4q, A\, . . . are defined via the formula 



.1 ; = 2 I — 



(355) 



for every j = 0, 1, 

In particular, the algorithm of Sections 15.21 15.31 15.41 with obvious minor modifications, is appli- 
cable to the task of evaluating PSWFs numerically. Moreover, the error analysis of such evaluation 
has been carried out in Theorem [3UJ and Section 15.51 

See Section [7] for numerical examples, involving the matrices similar to (|354[) . The results of 
some numerical experiments, where slightly modified versions of the algorithms of this paper (see 
Section [5]) are used to evaluate PSWFs, see, for example, [13], [15] . 



7 Numerical Results 

In this section, we illustrate the analysis of Section @] via several numerical experiments. All the 
calculations were implemented in FORTRAN (the Lahey 95 LINUX version), and were carried out 
in double precision. Occasionally, extended precision calculations were used to estimate the accuracy 
of double precision calculations. 

7.1 Experiment 1. 

In this experiment, we illustrate the performance of the algorithm on certain matrices. 
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c 
n 
A 


100 

250 
0.51665E+01 
0.37636E-39 


1,000 

2,100 
0.41168E+01 
0.29308E-41 


10,000 

20,215 
0.40117E+01 
0.23304E-42 


100,000 
200,500 
0.40012E+01 
0.37902E-43 


re/(Xi) 
rd(Yi) 


0.16471E-13 
0.10619E-13 


0.36528E-13 
0.58053E-13 


0.64311E-11 
0.67080E-12 


0.12414E-09 
0.94020E-11 


max(ois(X)) 
max(abs(l A )) 
max(rei(X)) 
max(reZ(Y)) 


0.23315E-14 
0.58981E-15 
0.27903E-12 
0.57672E-13 


0.12386E-14 
0.13891E-14 
0.26192E-10 
0.32331E-10 


0.27500E-12 
0.94943E-14 
0.77438E-07 
0.15776E-07 


0.27399E-11 
0.64567E-13 
0.43275E-03 
0.54285E-04 


avg(abs(X)) 
a.vg(abs(Y)) 
&vg{rel{X)) 
avg(rd(Y)) 


0.63819E-15 
0.74922E-16 
0.19247E-13 
0.53559E-14 


0.48374E-15 
0.61229E-15 
0.11704E-12 
0.15563E-12 


0.40466E-13 
0.43697E-14 
0.27750E-10 
0.59087E-11 


0.31112E-12 
0.25992E-13 
0.30753E-08 
0.39322E-09 



Table 1: Inverse Power (30 iterations) vs Principal Algorithm. 



Description. We proceed as follows. First, we choose, more or less arbitrarily, a real number 
c > 1. Then, we choose an integer number n ~ 2c + 10\/c, and construct a symmetric tridiagonal 
real n by n matrix A (see Definition [l}, whose diagonal entries A\, . . . ,A n are defined via (|292|) of 
Theorem [3D] with a = 2. In other words, 

f 

Aj = 2 + —, (356) 

for every j = 1, . . . , n. Then, we define the real number A via the formula 

4 • 40 

A = 4+ log(10), (357) 

7T 

and use A as an initial approximation to an eigenvalue of A for Inverse Power Method (see Sec- 
tion I3.4.ip . We use 30 iterations of Inverse Power Method to evaluate the eigenvalue A of A 
(A is closer to A than any other eigenvalue of A), and the corresponding unit length eigenvector 
Y = (Yi, . . . ,Y n ) 6 R", whose first coordinate is positive. We repeat the calculation in extended 
precision to obtain the vector Y ext . 

Next, we run the principal algorithm (see Sections 15.21 15.31 15. 4[) to evaluate the unit length 
A— eigenvector X = (Xi, . . . ,X n ) E K™ of A, whose first coordinate is positive (obviously, in exact 
arithmetics, we would have X = Y). We repeat the calculation in extended precision to obtain the 
vector X ext . 

We repeat the experiment four times, with c = 10 2 , 10 3 , 10 4 , 10 5 . First, we make the following 
observations. 

Observation 1. For every c = 10 2 , 10 3 , 10 4 , 10 5 , the eigenvalue A is evaluated to at least 15 
correct decimal digits (in double precision). In other words, the evaluated value of A is exact up to 
the machine precision. 

Observation 2. For every c = 10 2 , 10 3 , 10 4 , 10 5 , every coordinate of Y ext coincides with the 
corresponding coordinate of X ext in at least 16 decimal digits. Moreover, coordinate- wise verification 
confirms that both X ext and Y ext approximate the corresponding eigenvector of A to at least 16 
decimal digits. 

It follows from Observation 2 that any of X ext , Y ext can be used as a "reference eigenvector" , to 
evaluate the accuracy of each coordinate of X, Y. 
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Table structure. We display the results of the experiment in Table [T] This table has the 
following structure. The first column contains the description of the evaluated quantity, and the 
subsequent four columns correspond, respectively, to c = 10 2 , 10 3 , 10 4 , 10 5 . The first four rows 
contain c, the dimensionality n, the eigenvalue A and the first coordinate of the eigenvalue X\. The 
next two rows contain the relative errors of X\ and Y\ (see (|324[) ). The subsequent two rows contain 
the maximal absolute errors (see (|323p 'l among all the coordinates of X,Y, respectively, i.e. 



max(abs(X)) — max {abs(Xj)} , max(abs(Y)) = max {abs(Yj)} . 

l<j<n 



(358) 



The next two rows contain the maximal relative errors (see ([324])) among all the coordinates of X, Y, 
respectively, i.e. 

m&x(rel(X)) = max {rel(Xj)} , max(reZ(y)) = max {rel(Yj)} . (359) 

The subsequent two rows contain the average absolute errors among all the coordinates of X, Y, 
respectively, i.e. 



avg 



. n 1 n 

(abs(Xj) = - ^absiXj), avg(abs{Y)) = -^]a6s(Y, 



(360) 



The last two rows contain the average relative errors among all the coordinates of X, Y , respectively, 
i.e. 



avg 



1 n i n 

(rel(X)) = -J^reliXj), avg(rel(Y)) = - 5>ei(Yj). 



(361) 



Figures. Also, we display the results of the experiment in Figures [TJ [21 El HI In each of the 
figures, the x— axis corresponds to the indices of the eigenvector, i.e. 1 < j ' < n. Thus, each figure 
contains a plot of a certain function of the indices. 

The indices are divided into several groups (see also Definition [3]) . The color of the plot varies 
between different groups. The groups are defined via (|306[) . p07p . The integers k < p < m of (|307[) 
are marked on the x— axis of each plot. Also, each group is described by how the corresponding 
coordinates have been computed by the algorithm of Sections 15.21 15.31 15.41 The description of 
different groups of indices is described in Table [2] 



indices 


A - Aj 


based on 


described in 


color 


1 < j < k 


A - Aj > 2 


Theorem [7] 


Sectionl5.3.1l 


blue 


k < j < k + l 


2 > A - Aj 


Theorem [7] 


Section |5.3.2| 


black 


k + l < j <p 


2 > A - Aj > 


Theorem [TjJ 


Section 15.3.31 


red 


p < j < m — r 


> A - Aj > -2 


Theorem [Ml 


Section 15.3.61 


magenta 


m — r < j < m 


A - Aj > -2 


Theorem [7J 


Section 15.3.51 


black 


m < j < n 


-2 > A - Aj 


Theorem [7J 


Section |5.3.4| 


blue 



Table 2: Division of the indices. The integers l<k<k + l<p<m — r < m < n are those of 
(TO . ([Mil- 



Figures [U [21 [21 HI correspond, respectively, to c = 10 2 , 10 3 , 10 4 , 10 5 . Each of the four figures 
contains six plots. 

On Figures 1(a) 2(a) 3(a) 4(a)[ we plot the coordinates Xj of X, on the linear scale (left), and 
on the logarithmic scale (right). 
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On Figurcs [T(b)||2(b)||3(b)H4(b)[ wc plot the relative and absolute errors of X 3 (see f3"24]) . <IT2~5\i ), 
on the logarithmic scale: relative error (left), and absolute error (right). 

On Figures l(c)| 2(c) 3(c) 4(c)| we plot the relative and absolute errors of Yj (see Q324p . (I323[) ). 
on the logarithmic scale: relative error (left), and absolute error (right). 

Remark. On Figure [31 roughly every 5th point is plotted; On Figure [H roughly every 50th 
point is plotted (for the sake of visual enhancement). 

In addition, in Figure [5] we plot the following two functions of k + I < j < p (see Table [5]). First, 
we plot \(Xj\, where, a.j is that of Theorem [T71 fsee also Section l5.3.3[) . Also, we plot Xj/2, where Xj 
are the coordinates of the eigenvector. We recall that, due to (|148[) of Theorem[T7J Xj/2 = di(aj), for 
every j = . . . ,p. Figure[5]contains two plots, corresponding, respectively, to c = 100, n = 250 

and c = 1000, n = 2100. 

The following observations can be made from Table [l] Figures [l] [2J [3l HI [5l and additional 
experiments by the authors. 

Observation 1. For every c, the coordinates Xj of X behave (as a function of j) as predicted 
by Theorems [TT1 IT31 [T71 |2~41 (see also Definition [2]) . In other words, Xj grow fast up until j = k, 
i.e. in the region of growth; then, oscillate between j = k and j = to, i.e. in the oscillatory region; 
and, finally, decay fast (in absolute value) as j increases from to to n, i.e. in the region of decay (see 
Figures [I(a)| [2(a)] [3(a)] [4(a) I . 

Observation 2. Relative errors of both Xj and Yj exhibit fairly regular behavior in the regions 
of growth and decay, as expected (see Sections 15.3.11 [573 .41 15.51 in particular (J33HJ)) - 

Observation 3. As opposed to Observation 2, the relative errors of both Xj and Yj oscillate 
in the oscillatory region, varying by several orders of magnitude. This observation is not surprising 
(see Sections HTTP IS73761 l5~5l in particular (l340)) . (p42jn . 

Observation 4. While, for all four values of c, the first coordinate X\ is rather small (of order 
10 -40 ), it is evaluated fairly accurately by both the principal algorithm and Inverse Power Method 
(see Section l5~5j) . 

Observation 5. According to the combination of Theorem 1301 (|343[) in Section [5751 and (|356p . 
we expect to lose roughly 2 decimal digits in the evaluation of X\, as we multiply c by the factor 
of 10. This prediction is confirmed by the data in Table Q] (the fifth row): as c changes from 10 3 to 
10 4 to 10 5 , the relative error increases from 0.4E-13 to 0.6E-11 to 0.1E-9. Still, even for as large c 
as 10 5 and as small X\ as w 10~ 40 , the latter is evaluated to 10 decimal digits. 

Observation 6. While the coordinates Xj oscillate rather rapidly in the oscillatory regime, the 
absolute value of the complex parameter (Xj (see Theorem fT7[) changes fairly slowly (see Figure [5j . 
This observation provides an intuitive justification for the steps of the principal algorithm, described 
in Sections [5.3.31 [573751 Also, it agrees with the error analysis of Section [5751 

Observation 7. The coordinates of the vector Y are usually evaluated somewhat more accu- 
rately than those of X (by roughly one decimal digit). In other words, 30 iterations of Inverse Power 
Method, while being more expensive by a constant factor in terms of operations, slightly outperform 
the principal algorithm of this paper. At the bottom line, however, the two algorithms are rather 
similar in terms of both computational cost and accuracy. 



7.2 Experiment 2. 

In this experiment, we illustrate the performance of the algorithm on the matrices similar to those 
of Experiment 1. 

Description. We proceed as follows. We choose the real number c = 1000 and integer n = 2100. 
Then, we construct the symmetric tridiagonal n by n matrix A, whose non-zero off-diagonal entries 
are all one, and whose diagonal entries A\, . . . , A n are defined via (|356|) in Section [77T1 

We use Inverse Power Method (see Section 13.4.11) to evaluate several eigenvalues of A, and, for 
each eigenvalue A, we evaluate the corresponding unit length eigenvector Y — Y(X) — (Yi, . . . , Y n ) e 
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M. n , whose first coordinate is positive (using 30 iterations of Inverse Power Method for each A). 
For each such eigenvalue A, we repeat the calculation in extended precision to obtain the vectors 

Next, for each such A, we run the principal algorithm (see Sections 15.21 15.31 15.41) to evaluate the 
unit length A— eigenvector X = X(X) = (X\, . . . ,X n ) € R™ of A, whose first coordinate is positive 
(obviously, in exact arithmetics, we would have X = Y). For each A, we repeat the calculation in 
extended precision to obtain the vectors X ext — X ext (X). 

First, we make the following observations. 

Observation 1. In double precision, all the evaluated eigenvalues A are correct to tat least 15 
decimal digits. In other words, these eigenvalues are exact up to the machine precision. 

Observation 2. For every A, every coordinate of Y ext (X) coincides with the corresponding 
coordinate of X ext (X) in at least 16 decimal digits. Moreover, coordinate- wise verification confirms 
that, for every A, both X ext (X) and Y ext (X) approximate the corresponding eigenvector of A to at 
least 16 decimal digits. 

It follows from Observation 2 that, for every A, any of X ext (X), Y ext (X) can be used as a "reference 
eigenvector", to evaluate the accuracy of each coordinate of -^(A), Y(X). 



X 




rel(Xi) 


max(reZ(X)) 


avg(reZpO) 





40351E- 


-01 





10809E- 


13 





58646E- 


12 





29176E-09 





13114E- 


11 





40471E- 


-01 





99452E- 


18 





66326E- 


12 





28637E-09 





17473E- 


11 





40595E- 


-01 





63720E- 


22 





86132E- 


12 





99517E-09 





23739E- 


11 





40705E+01 





12641E- 


25 





28375E- 


11 





30593E-08 





81554E- 


11 





40836E- 


-01 





46754E- 


30 





10490E- 


12 





13829E-09 





63507E- 


12 





40932E- 


-01 





27309E- 


33 





70170E- 


12 





13838E-09 





16682E- 


11 





41069E+01 





66341E- 


38 





35095E- 


12 





59918E-09 





14458E- 


11 





41168E+01 





29308E- 


41 





36528E- 


13 





26192E-10 





11704E- 


12 





41289E- 


-01 





24013E- 


45 





42306E- 


12 





29842E-08 





27627E- 


11 





41392E+01 





84484E- 


49 





87444E- 


13 





15693E-09 





60358E- 


12 





41537E- 


-01 





10509E- 


53 





70798E- 


12 





25694E-09 





12596E- 


11 





41643E- 


-01 





29507E- 


57 





17023E- 


13 





62182E-09 





54695E- 


12 





41728E- 


-01 





39899E- 


60 





29785E- 


11 





76945E-10 





39780E- 


11 


0.42665E- 


-01 


0.13675E-91 


0.34374E-13 


0.20497E-09 


0.60583E-12 



Table 3: Principal Algorithm. Parameters: c — 1000, n — 2100. 

Table structure. The results of the experiment are displayed in Tables [31 0J These tables 
have the following structures. The first column contains some eigenvalues A of the matrix A. The 
second column contains the first coordinate of the eigenvector X (Table [3]) or Y (Tabled]). The 
third column contains the relative error of this coordinate (see (1324j) '). The fourth column contains 
the maximal relative error of X, Y, respectively (see (|359[) in Section ITTTj) . The last column contains 
the average relative error of X, Y, respectively (see (|36ip in Section l7T2"j) . 

The following observations can be made from Tables [21 [H 

Observation 1. Usually, the vector Y(X) is more accurate than the vector X(X) by roughly one 
decimal digit. 

Observation 2. The accuracy of X(X) is roughly the same throughout all of Tableau in terms of 
relative error. In particular, the first coordinate X\ of the eigenvector X is evaluated by the principal 
algorithm with roughly the same accuracy, for every A, independent of the order of magnitude of 
X\. In other words, even though X\ varies between ss 10~ 14 (first row) and 10 -92 (last row), it is 
always evaluated to roughly 12 decimal digits (see also Section [533J • 
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rel(Yx) 


max(reZ(Y")) 


avg(rel(Y)) 


403 c >1 F+01 


U 


1 nonoTT 
lUoUyrj- 




U.z4yoyiir 


lo 


U 


ioZoZrj- 


no 

uy 


U 






o 4n47i F+01 


n 
u 




1 R 




±o 


n 

u 




1U 


n 
u 




1 9 


0.40595E+01 





63720E- 


22 


0.54421E- 


13 





15798E- 


09 





26474E- 


12 


0.40705E+01 





12641E- 


25 


0.49833E- 


13 





12447E- 


09 





26264E- 


12 


0.40836E+01 





46754E- 


30 


0.47955E- 


13 





30739E- 


10 





16290E- 


12 


0.40932E+01 





27309E- 


33 


0.63108E- 


13 





35663E- 


10 





15882E- 


12 


0.41069E+01 





66341E- 


38 


0.40131E- 


13 





60993E- 


10 





19122E- 


12 


0.41168E+01 





29308E- 


41 


0.58053E- 


13 





32331E- 


10 





15563E- 


12 


0.41289E+01 





24013E- 


45 


0.61224E- 


13 





53569E- 


09 





40129E- 


12 


0.41392E+01 





84484E- 


49 


0.48218E- 


13 





50575E- 


10 





21127E- 


12 


0.41537E+01 





10509E- 


53 


0.64663E- 


13 





29615E- 


10 





13552E- 


12 


0.41643E+01 





29507E- 


57 


0.49151E- 


13 





27351E- 


09 





25378E- 


12 


0.41728E+01 





39899E- 


60 


0.48311E- 


13 





83383E- 


11 





11069E- 


12 


0.42665E+01 


0.13675E-91 


0.46828E-13 


0.33287E-10 


0.74736E-13 



Table 4: Inverse Power (30 iterations) . Parameters: c = 1000, n = 2100. 



Observation 3. The accuracy of Y(X) is roughly the same throughout all of TablelH in terms of 
relative error. In particular, the first coordinate Yx of the eigenvector Y is evaluated by the principal 
algorithm with roughly the same accuracy, for every A (to roughly 13 decimal digits), independent 
of the order of magnitude of Yx . 

7.3 Experiment 3. 

In this experiment, we illustrate the numerical algorithms of Section [5] via evaluation of Bessel 
functions (see Sections 13.21 13. 4. 31 16. ip . 

Description. We choose, more or less arbitrarily, the real number c > and the integer n > 0. 
Then, we choose the integer N > max {n,c}. We construct the symmetric tridiagonal 2N + 1 by 
2N + 1 matrix A, whose non-zero off-diagonal entries are all one (see (|6"0)l ). and whose diagonal 
entries Ax, ... , A2N+X are defined via (|6"Tj) of Remark [5] in Section f3. 4.31 We define the real number 
A via fl6"3"]) (A is an eigenvalue of A, due to Remark 

We use the principal algorithm (see Sections 15.21 15.31 15. 4j) to evaluate the A— eigenvector X E 
K 2JV+1 of A whose first coordinate is positive. We repeat the computation in extended precision, to 
obtain X ext £ M. 2N+1 . Then, we run 30 iterations of Inverse Power Method (see Section r3.4.1j) to 
evaluate the unit length A— eigenvector Y € M 2Ar+1 of A whose first coordinate is positive (obviously, 
in exact arithmetics, X and Y would be the same vector). We repeat the computation in extended 
precision to obtain Y ext £ R 2N+1 . 

Based on Remark[9l we evaluate Jo(c), . . . , J n (c) from X (or Y) as follows. Suppose, for example, 
that 

X = (Xjf, Xjsr-Xi • • • j ^i>^0) X-x, • • • , X_(jv-i),X_jv). (362) 

We evaluate Jo(c) as X , Ji(c) as Xx, up to J„(c) as X n . 

We repeat the experiments for all the choices of parameters c, n, N, listed in Table [5] 
First, we make the following observations. 

Observation 1. For c = 100, we ran the experiment in extended precision with N = 215 and 
N = 235, to obtain X ext G R 2JV+1 and Y ext G WL 2N+1 , for each N. These four vectors coincide in the 
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c 


n 


N 


100 


200 


215 


100 


200 


235 


1,000 


1,200 


1,250 


1,000 


1,200 


1,270 


10,000 


10,490 


10,550 


10,000 


10,490 


10,570 


100,000 


101,000 


101,150 


100,000 


101,000 


101,200 



Tabic 5: Choice of parameters for Experiment 3. 

middle 2n + 1 = 401 coordinates up to at least 15 digits. In other words, using the notation (|362[) . 
these four vectors coincide in the coordinates with indices 0, ±1, . . . , ±n up to at least 15 digits. 

Observation 2. Similar to Observation 1, for each c = 10 3 , 10 4 , 10 5 and each corresponding 
pair of N from Table [5l the four vectors coincide in the middle 2n + 1 coordinates up to at least 15 
digits (in extended precision). 

We conclude from Observations 1, 2 that, for each pair of c,n from Table for either of the 
two corresponding values of N, both X ext and Y ext can be used as "reference" vectors to evaluate 
the accuracy of the evaluation of Jo(c), . . . , J n (c). In particular, in extended precision, each of 
Jo(c), . . . , J„(c) is evaluated to at least 15 decimal digits, by either method. 



c 


N 


k 


J fe (c) 


\X k /J k (c) - 1| 


\Y k /J k (c) - 1| 


100 


215 





0.19986E-01 


0.48607E-14 


0.14408E-13 






1 


-.77145E-01 


0.69438E-13 


0.17989E-14 






2 


-.21529E-01 


0.83800E-14 


0.13053E-13 






200 


0.20594E-40 


0.48517E-13 


0.29705E-14 


1,000 


1,250 





0.24787E-01 


0.27994E-14 


0.83983E-15 






1 


0.47283E-02 


0.10944E-11 


0.84382E-13 






2 


-.24777E-01 


0.32206E-14 


0.00000E+00 






1,200 


0.83509E-38 


0.11268E-12 


0.53135E-14 


10,000 


10,550 





-.70962E-02 


0.17552E-12 


0.25069E-12 






1 


0.36475E-02 


0.13038E-10 


0.93087E-12 






2 


0.70969E-02 


0.17648E-12 


0.25042E-12 






10,490 


0.35152E-46 


0.14137E-11 


0.16278E-12 


100,000 


101,150 





-.17192E-02 


0.76153E-11 


0.88710E-11 






1 


0.18468E-02 


0.48246E-10 


0.76843E-11 






2 


0.17192E-02 


0.76138E-11 


0.88708E-11 






101,000 


0.39770E-43 


0.28153E-10 


0.11461E-11 



Table 6: Evaluation of J k (c). Corresponds to Experiment 3. 

Next, we run the experiment in double precision, with c = 10 2 , 10 3 , 10 4 , 10 5 , and N = 215, 1250, 10550, 101150, 
respectively. 

Table structure. The results of this experiment are displayed in Table [5] This table has the 
following structure. The first two columns contain the parameters c and N . The third column 
contains the integer k between and n, where n < N is that of Table [5] The fourth column contains 
J k {c) (evaluated in extended precision). The fifth column contains the relative error to which X k 
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approximates Jfc(c). The last column contains the relative error to which Yk approximates Jfc(c). 

We make the following observations from Tableland some additional experiments. 

Observation 1. For every c, the relative accuracy of X n is similar to that of Xq, Xi, X2. This 
is despite the fact that \Xk\ ~ 10~ 2 for k = 0, 1, 2, while \X n \ ps 10~ 40 . This phenomenon, however, 
confirms the analysis of Section l5~5l 

Observation 2. According to (|343[) of Section I5~51 (see also Section T3.4.3p . when we evaluate 
J„(c) via computing X n , we expect to lose roughly 4/3 decimal digits every time we increase c by 
the factor of 10. This prediction is confirmed by the data in Table [S] (see rows 4,8,12,16 in the fifth 
column) . 

Observation 3. The coordinate of Y, evaluated via 30 iterations of Inverse Power Method, are 
usually more accurate than those of X (by roughly one decimal digits). However, sometimes the 
opposite is true (see rows 1,5,9,13). 
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x., 0=100 iog 10 (|x.|), 0=100 




k=108 p=176 m=227 k=108 p=176 m=227 



(a) coordinates: linear and logarithmic scales 
log 10 (rel(X.)), c=100 log 10 (abs(X.)), c=100 




k=108 p=176 m=227 k=108 p=176 m=227 



(b) principal algorithm: relative and absolute errors 
log 1() (rel(Y.)), c=100 log 1f) (abs(Y)), c=100 




k=108 p=176 m=227 k=108 p=176 m=227 



(c) inverse power: relative and absolute errors 

Figure 1: The coordinates of X (principal algorithm) and Y (30 iterations of Inverse Power). 
Parameters: c = 100, n = 250, A = 0.51665E+01. 
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X., c=1000 



0.05 



-0.05 




log 10 (|X.|),c=1000 




\ 



k=341 



p=1453 m=2028 



k=341 



p=1453 m=2028 



(a) coordinates: linear and logarithmic scales 
log 10 (rel(X.)), c=1000 log 10 (abs(X.)), c=1000 





mm 



\ 



k=341 



p=1453 m=2028 



k=341 



p=1453 m=2028 



(b) principal algorithm: relative and absolute errors 



log 10 (rel(Y.)),c=1000 



log 10 (abs(Y.)), c=1000 





k=341 



p=1453 m=2028 



(c) inverse power: relative and absolute errors 



Figure 2: The coordinates of X (principal algorithm) and Y (30 iterations of Inverse Power). 
Parameters: c = 1000, n = 2100, A = 0.41168E+01. 



59 



X., 0=1 0,000 




log 10 (|X.|),c=10,000 



p=14182 m=20029 




p=14182 m=20029 



(a) coordinates: linear and logarithmic scales 



log 1Q (rel(X.)), c=1 0,000 




k=1082 



p=14182 m=20029 



log 10 (abs(X.)), c=1 0,000 




k=1082 



p=14182 m=20029 



(b) principal algorithm: relative and absolute errors 



log 1Q (rel(Y.)), 0=1 0,000 




-15L _ 

k=1082 



p=14182 m=20029 



log 10 (abs(Y.)), c=1 0,000 




k=1082 



p=14182 m=20029 



(c) inverse power: relative and absolute errors 



Figure 3: The coordinates of X (principal algorithm) and Y (30 iterations of Inverse Power). 
Parameters: c = 10, 000, n = 20, 215, A = 0.40117E+01. 
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X., c=1 00,000 



log 10 (|X.|),c=100,000 



0.02 


1 









0.01 


1 




-10 



















-20 














t 




-30 












0.01 




" 4 ° 













k=3423 p=141461 m=200029 k=3423 p=141461 m=200029 



(a) coordinates: linear and logarithmic scales 
log 10 (rel(X.)), c=100,000 log 10 (abs(X.)), c=100,000 




k=3423 p=141461 m=200029 k=3423 p=141461 m=200029 



(b) principal algorithm: relative and absolute errors 



log 10 (rel(Y.)), c=1 00,000 



-10 
-12 
-14L 




log 10 (abs(Y.)), c=1 00,000 



-15 
-20 

-30 

-40 

-50 



k=3423 



p=141461 m=200029 k=3423 



p=141461 m=200029 



(c) inverse power: relative and absolute errors 



Figure 4: The coordinates of X (principal algorithm) and Y (30 iterations of Inverse Power). 
Parameters: c = 100, 000, n = 200, 500, A = 0.40012E+01. 
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(a) c = 100, n = 250, A = 0.51665E+01 
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(b) c = 1, 000, n = 2, 100, A = 0.41168E+01 
Figure 5: The coordinates Xj vs. \ctj\. Corresponds to Experiment 1. 
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